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ABSTRACT 

We have combined the semi-analytic galaxy formation model of Guo et al. (2011) with the 
particle-tagging technique of Cooper et al. (2010) to predict galaxy surface brightness pro- 
files in a representative sample of ~ 1900 massive dark matter haloes (10 12 -10 14 Mq) from 
the Millennium II ACDM N-body simulation. Here we present our method and basic results 
focusing on the outer regions of galaxies consisting of stars accreted in mergers. These simu- 
lations cover scales from the stellar haloes of Milky Way-like galaxies to the 'cD envelopes' of 
groups and clusters, and resolve low surface brightness substructure such as tidal streams. We 
find that the surface density of accreted stellar mass around the central galaxies of dark matter 
haloes is well described by a Sersic profile, the radial scale and amplitude of which vary sys- 
tematically with halo mass (M2oo)- The total stellar mass surface density profile breaks at the 
radius where accreted stars start to dominate over stars formed in the galaxy itself. This break 
disappears with increasing M200 because accreted stars contribute more of the total mass of 
galaxies, and is less distinct when the same galaxies are averaged in bins of stellar mass, be- 
cause of scatter in the relation between M* and M 2 oo- To test our model we have derived 
average stellar mass surface density profiles for massive galaxies at z « 0.08 by stacking 
SDSS images. Our model agrees well with these stacked profiles and with other data from 
the literature, but can be more rigorously tested by future surveys that extend the analysis of 
the outer structure of galaxies to fainter isophotes. We conclude that it is likely that the struc- 
ture of the spheroidal components of galaxies is largely determined by collisionless merging 
during their hierarchical assembly. 

Key words: galaxies: structure; galaxies: elliptical and lenticular, cD; galaxies: bulges; galax- 
ies: evolution 



1 INTRODUCTION 

Hierarchical clustering leads to the coalescence of dark matter 
haloes. Galaxies formed 'in situ' by dissipative collapse in the cores 
of these haloes (White & Rees 1978) accrete additional stars from 
the debris of their hierarchical progenitors. The aim of this paper 
is to predict how the surface density profiles of galaxies reflect 
changes in the balance between in situ star formation and stellar 
accretion during their hierarchical growth over the lifetime of the 
universe. 

The idea that observations of accreted stars could be used 
to test theories of galaxy evolution has its roots in the study of 
the stellar halo and globular clusters of the Milky Way and M31 
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(Baade 1944; Eggen, Lynden-Bell & Sandage 1962; Searle & Zinn 
1978). The recent discovery of cold stellar streams in these haloes, 
some with identifiable progenitors, has provided direct evidence 
that they grow at least partly through the tidal disruption of com- 
panion galaxies (e.g. Ibata, Gilmore & Irwin 1995; Belokurov et al. 
2006; Niederste-Ostholt et al. 2010; McConnachie et al. 2009). 
Stellar haloes and streams appear to be a generic feature of late- 
type galaxies (Zibetti, White & Brinkmann 2004; Richardson et al. 
2009; Martmez-Delgado et al. 2010a; Bailin et al. 2011; Radburn- 
Smith et al. 2011). Shell-like structures have been detected around 
both early and late type galaxies (Malin & Carter 1983; Schweizer 
1980; Schweizer et al. 1990; Schweizer & Seitzer 1992; Tal et al. 
2009; Martmez-Delgado et al. 2010a) and can also be readily ex- 
plained as the result of galactic accretion in a cold dark matter 
(CDM) universe (e.g Cooper et al. 201 1). Galaxies at the centres of 
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massive clusters are often surrounded by extended envelopes of dif- 
fuse 'intracluster light' (ICL) (Matthews, Morgan & Schmidt 1964; 
Oemler 1976; Thuan & Romanishin 1981; Schombert 1988; Gra- 
ham et al. 1996; Lin & Mohr 2004; Gonzalez, Zabludoff, & Zarit- 
sky 2005; Mihos et al. 2005; Krick, Bernstein, & Pimbblet 2006; 
Donzelli, Muriel, & Madrid 2011) which is also thought to orig- 
inate from the stripping and disruption of satellite galaxies (Gal- 
lagher & Ostriker 1972). 

Semi-analytic models of galaxy formation aim to quantify the 
importance of the accretion of stars and gas in different types of 
galaxy (White & Frenk 1991; Cole 1991; Kauffmann, White & 
Guiderdoni 1993; Somerville & Primack 1999; Cole et al. 2000; 
Baugh et al. 2005; De Lucia & Blaizot 2007; Benson & Bower 
2010; Guo et al. 201 1). In particular, they predict how the mass of 
stars accreted by a galaxy depends on the mass and assembly time 
of its dark matter halo, as well as the number of progenitors of the 
halo and their individual star formation histories. Such models have 
shown that stellar accretion and merging can supply the majority of 
stars in the most massive galaxies at the present day, but make only 
a minor contribution in galaxies with masses less than that of the 
Milky Way, where most stars form from gas that cools directly from 
the halo (Baugh, Cole & Frenk 1996; Kauffmann 1996; De Lucia 
et al. 2006; Naab et al. 2007; Purcell, Bullock, & Zentner 2007; 
Guo & White 2008; Parry, Eke, & Frenk 2009). 

Accretion and merger events may affect the colour, size and 
morphology of a galaxy, even if they make a limited contribution to 
its mass. These effects are thought to be relevant to the dichotomy 
between early and late-type morphologies in the Hubble sequence 
(e.g. Toomre 1977; ; Frenk et al. 1985; Cowie et al. 1994; Zepf 
1997; Kauffmann & Chariot 1998; Cole et al. 2000; Bell et al. 2004; 
Sales et al. 2012) and many well-known scaling relations between 
observable properties of massive galaxies (Faber & Jackson 1976; 
Kormendy 1977; Djorgovski & Davis 1987; Peletier et al. 1990; 
Bender, Burstein & Faber 1992; Kauffmann et al. 2003b; Bernardi 
et al. 2003), including correlations between the luminosity of mas- 
sive elliptical galaxies and the extent and shape of their projected 
surface brightness profiles (Kormendy 1977; Binggeli, Sandage & 
Tarenghi 1984; Schombert 1986; Graham et al. 1996; Graham & 
Guzman 2003; Kormendy et al. 2009). 

The effects of stellar accretion on galactic structures depend 
not only on the population of infalling galaxies and the rate at which 
haloes coalesce, but also on the gravitational dynamics of the accre- 
tion process. This means that galactic accretion cannot be studied in 
isolation from galaxy formation. The way in which accreted stars 
are deposited in the outer regions of galaxies and the consequent 
change in observables such as half-light radius and stellar mass 
has been considered extensively in recent literature (Daddi et al. 
2005; Trujillo et al. 2006; van Dokkum et al. 2010). Simulations 
focussing on this issue by e.g. Naab, Khochfar & Burkert (2006); 
Naab et al. (2007); Naab, Johansson, & Ostriker (2009) and Oser 
et al. (2010) have highlighted the importance of N-body dynam- 
ics when making specific predictions for the evolution of galaxy 
sizes and velocity dispersions in CDM (compare Gonzalez et al. 
2009; Guo et al. 201 1; Shankar et al. 2013). Notably, these simula- 
tions suggest that high mass ratio mergers contribute significantly 
to the structure of massive elliptical galaxies (Hilz, Naab & Ostriker 
2012a; Hilz et al. 2012b). 

Galaxy surface brightness or surface mass density profiles 
can be described by Sersic (1968) functions (Caon, Capaccioli & 
D'Onofrio 1993; Graham & Guzman 2003; Graham 201 1, and ref- 



erences therein) of the form 

E(i?) = S 50 e X p{-b n [(R/R 50 ) 1/n - 1]}, (1) 

where Rso is the 'effective' radius which encloses half the total 
mass or light in projection, E50 = £(i?5o) is the surface density at 
7?50, which sets the amplitude of the profile, and n is the Sersic in- 
dex, which sets the concentration of the profile. Profiles with index 
n — 1 are exponential and thus a good fit to late type (disk) galax- 
ies, dwarf elliptical galaxies (Binggeli & Cameron 1991; Graham 
& Guzman 2003) and some (not all) galactic bulges (Andredakis, 
Peletier & Balcells 1995). Galaxies with early type morphology 
typically have Sersic index n > 2.5, with a characteristic value of 
n = 4 for massive ellipticals. In general n increases with luminos- 
ity such that some of the central galaxies of massive clusters have 
approximately power-law profiles (n > 10) (e.g. Graham et al. 
1996). Collisionless numerical simulations predict that the relaxed 
remnants of binary galaxy mergers have n ~ 4 profiles (Negro- 
ponte & White 1983; Aguilar & White 1986; Barnes 1988). These 
results (and an emphasis on n = 4 profiles in observational stud- 
ies) have given rise to a frequently encountered 'working model' 
for elliptical galaxies. In this model, deviations of surface bright- 
ness from the de Vaucouleurs (1948) (n = 4 Sersic) profile are 
interpreted as evidence for additional (or even 'missing') structural 
components, separate from the central galaxy itself. 

In this paper we use an extension of the Guo et al. (2011, 
hereafter Gl 1) semi-analytic galaxy formation model to predict the 
spatial distribution of all stars accreted on to present-day massive 
galaxies. Like many recent numerical studies of the size evolution 
of massive quiescent galaxies (Meza et al. 2003; Naab et al. 2009; 
Oser et al. 2010) and Milky Way-like galaxies (Abadi, Navarro & 
Steinmetz 2006; Cooper et al. 2010; Font et al. 2011; Sales et al. 
2012) we emphasize the difference between in situ star formation 
and galactic accretion in our analysis. The semi-analytic model al- 
lows us to compute the full in situ star formation histories of our 
galaxies and all their hierarchical progenitors. In addition, we use 
an N-body method to predict how each accreted population evolves 
in all six dimensions of phase space. We demonstrate that correla- 
tions between Sersic parameters and galaxy luminosity can be ex- 
plained by changes in the relative proportion of stars formed in situ 
and stars accreted as a function of halo mass. Because our model 
is applied to the Millennium II simulations, which contain ~ 2000 
haloes in the mass range 10 12 < M200 < 1O 14 M , we can make 
a statistical comparison to observational data on the total amount 
of accreted light around galaxies and its distribution, using simu- 
lations that match other constraints such as the galaxy stellar mass 
function. In this paper, we only consider galaxies at the centres of 
virialised dark haloes at z = 0, although our model treats the entire 
galaxy population across cosmic time and can thus make predic- 
tions for galaxies at high redshift and for galaxies that are satellites 
at the present day. We focus on the spatial distribution of accreted 
stars, but our model also predicts their kinematic properties and 
chemical abundances. 

A statistical study of this sort is motivated by the availability 
of moderately deep wide-field imaging from surveys such as SDSS 
Stripe 82 and PanSTARRS, which will enable us to determine the 
extent to which results obtained from the Milky Way and M3 1 stel- 
lar haloes are applicable to the galaxy population as a whole. Previ- 
ous studies of galaxy structure in large surveys have focused on re- 
gions of high surface brightness (e.g. Shen et al. 2003), because the 
outskirts of galaxies are usually much fainter than the sky, even in 
the case of the envelopes around brightest cluster galaxies (BCGs). 
Stacking (e.g. Zibetti et al. 2005; Tal & van Dokkum 2011) is a 
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promising technique for studying the average properties of these 
regions. We have therefore carried out our own stacking analysis 
using imaging data from SDSS, for comparison to our models. 

We proceed as follows. In Section [2] we describe how we se- 
lect a sample of massive central galaxies from the Millennium II 
simulation. We also summarize how our particle-tagging method 
works. Section[3]shows examples of the stellar haloes of individual 
galaxies in our model. In Section|4]we present our main statistical 
results in the form of average stellar mass surface density profiles 
for haloes and galaxies of different mass, highlighting differences 
between in situ and accreted stars. We compare with observations 
from the literature in Section|5]and from our own stacking analysis 
in Section [6] Section [7] interprets trends in surface brightness pro- 
file shape by studying the origin of accreted stars in our simulation. 
We summarise our findings in Section [5] Appendix [A] discusses 
the numerical convergence of our method and differences with our 
previous work on particle tagging. Appendix [B] describes a simple 
numerical method that can be used to approximate the distributions 
of in situ stars at z = in our model. Appendix [C] describes the 
technical details of our SDSS stacking analysis. 



2 SIMULATIONS 

2.1 Semi-analytic model 

Millennium II (Boylan-Kolchin et al. 2009) is a collisionless N- 
body simulation of ACDM structure formation in a comoving vol- 
ume of 10 6 h~ 3 Mpc 3 with a flat ACDM cosmology, ft m = 0.25, 
Qa ~ 0.75 and Hubble parameter h = 0.73. The particle mass res- 
olution is 6.89 x 10 6 h' 1 M Q . The semi-analytic galaxy model of 
Gl 1 is based on halo merger trees derived from Millennium II and 
its parameters are tuned to fit the SDSS stellar mass function of ( ; 
shown as a grey dashed line in Fig.[TJ. Figure 7 of Gl 1 shows that 
when same model is applied to the Millennium simulation (which 
has a larger volume), it overpredicts the number of galaxies with 
log 10 M*/Mq > 11.5. Gil suggest that this discrepancy is due 
to sample variance and ~ 1 mag luminosity uncertainties for the 
most luminous galaxies in SDSS. 

From the results of Gl 1 we select 1872 central galaxies more 
massive than Af* = 5 x 10 10 Mq at z = 0. The stellar mass func- 
tion of these galaxies is shown by the black histogram in Fig.[TJ We 
do not select all the galaxies in the simulation above our threshold 
stellar mass (in particular, we exclude satellite galaxies, the central 
galaxy of the most massive cluster and a small number of central 
galaxies erroneously classified as its satellites). The grey histogram 
shows the mass function of all galaxies in the simulation to demon- 
strate that this selection does not bias our sample. 

The Gil model uses a combination of two ID axisymmetric 
density profiles to represent the distribution of stars inside galax- 
ies (an exponential disc and a Jaffe-model bulgQ and quantifies 
galaxy morphology using the ratio of bulge mass to total stellar 
mass, B/T. The red histogram in Fig. [T] shows the mass function 
of galaxies with B/T ^ 0.9 ('ellipticals') and the blue histogram 
the mass function of galaxies with B/T < 0.9. We plot the ratio 



1 Although Gil also track the total mass in a diffuse stellar halo component 
(as described in appendix A) they do not specify its density profile. The 
mass in this component is only a significant fraction of the total galaxy 
stellar mass above log 10 M* /Mq ~ 11.2. Horizontal bars in Fig.[T]show 
how including the Gil stellar halo component in the central galaxy stellar 
mass affects the stellar mass function at z = 0. 
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Figure 1. The solid black histogram shows the mass function of our primary 
sample of central galaxies from the Gil model. The solid grey histogram is 
the mass function of all galaxies in the Millennium II simulation, including 
satellite galaxies and all galaxies associated with the merger tree of the most 
massive cluster which we have excluded from our analysis. The dashed grey 
line shows the Li & White (2009) SDSS mass function normalized to the 
Millennium II volume, derived as in Appendix A of Guo et al. (2010) and 
assuming the same Chabrier IMF as our model. Red and blue histograms 
split the mass function of our primary sample into two components based 
on bulge-to-total mass ratio, as shown in the legend. Short black horizontal 
lines (between the black and grey histograms) show the effect of includ- 
ing the Gil predictions for diffuse stellar halo mass in each galaxy when 
computing the mass function of our sample. Considering only our sample 
galaxies, the lower two panels show (top) the fraction, / e n, of galaxies in 
each mass bin with B/T > 0.9, and (bottom) the fraction of the combined 
bulge and disc stellar mass in each bin that was formed in situ, split by B/T 
as in the main panel. 



of these mass functions in the middle panel. The fraction of 'el- 
liptical' galaxies increases from 50 to 100 per cent in the interval 
11.0 < log 10 M 4 /Mq < 11.3. This is in good agreement with ob- 
servations (e.g. Conselice 2006, see figure 4 of Gil) and a similar 
transition scale has been found in other semi-analytic models (De 
Lucia et al. 2011). These models also predict differences between 
early and late-type K-band luminosity functions that are qualita- 
tively similar to those observed (Benson & Devereux 2010). 

We define 'in situ' stars as those that are still gravitationally 
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bound to the dark matter halo in which they formed^] Gil predict 
that more massive galaxies form less of their total stellar mass in 
situ. The lowest panel of Fig. [T] shows that the fraction of stars in 
each mass bin formed in situ decreases from almost (but of course 
not exactly) 100 per cent at log 10 M*/Mq = 10.7 to 50 per cent 
at log 10 M*/Mq = 11.1. The mass fraction of in situ stars in sys- 
tems with B/T ^ 0.9 is similar to that of other galaxies regardless 
of stellar mass. This implies that in the Gil model, the relative 
contribution of in situ star formation depends primarily on stellar 
mass and not morpholog)]^] For the most massive galaxies plotted 
Fig. [T| the fraction of stars formed in situ is ~19 per cent. Thus 
the model makes a clear prediction that accreted stellar populations 
will dominate the structure of the most massive galaxies. 

2.2 Particle tagging 

We use a technique we call particle tagging to predict the stellar 
population mix and spatial distribution of stars in galaxies, based 
on the merger trees and star formation histories of the Gil model. 
This technique uses additional information from the underlying N- 
body simulation in order to predict more observables than stan- 
dard semi-analytic models, without running a new simulation. Ear- 
lier studies using related techniques include Bullock, Kravtsov & 
Weinberg (2001), Penarrubia, Navarro & McConnachie (2008) and 
Bullock & Johnston (2005), although these were not coupled to the 
predictions of semi-analytic models. We give a brief summary of 
our method below; for more detail see Cooper et al. (2010, hereafter 
C10). A discussion of minor differences between our implementa- 
tion and that of C10 and a test of convergence with their results can 
be found in Appendix |A| 

The particle tagging technique associates ('tags') sets of dark 
matter particles in an N-body simulation (here Millennium II) with 
stellar populations of a single metallicity and age. The tagged par- 
ticles can be used to track the evolution of their associated popula- 
tion in phase space, from the time when the stars form to the present 
day (z = 0). Our definition of a stellar population comprises all the 
stars formed in a single galaxy between two consecutive snapshots 
of the Gl 1 model. An isolated galaxy that forms stars at a constant 
rate for a Hubble time will produce a number of these populations 
equal to the number of simulation snapshots. All model galaxies at 
z = are a superposition of many such populations, because they 
accrete populations formed in their hierarchical progenitors as well 
as forming their own stars in situ. 

For every population, particles are selected according to a tag- 
ging criterion (described below). An equal fraction of the total mass 
of the population is given to each particle thus selected. Every new 
population tags a new set of particles, selected from the correspond- 
ing dark matter halo at the snapshot immediately after the popula- 
tion forms. This means that a DM particle can be tagged more than 
once, if it meets the tagging criterion for two or more populations 
(by construction, this can only happen at different snapshots). In 
such cases, each tag is tracked separately. The corollary of this is 
that each tagged particle carries its own unique star formation his- 
tory, with the time resolution of the Millennium II snapshots. 

2 The definition of a dark matter halo and the construction of halo merger 
trees follow De Lucia et al. (2006). 

3 Because galaxies with B/T > 0.9 were formed by low mass ratio merg- 
ers, one might expect that their in situ fraction would be lower than that of 
galaxies at the same mass without lower B/T > 0.9. However, in the Gil 
model, merger-induced starbursts and disk instabilities can increase the in 
situ mass while also increasing B/T. 
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Figure 2. Dots show in situ surface density profiles in two Milky Way-like 
haloes at z = from Gil, predicted by our particle tagging model with 
/mb = 1% (blue) and 10% (red). Upper and lower panels respectively cor- 
respond to galaxies with M 2 qo = (12.1, 12.3), M* = (10.8, 10.9) and 
NFW concentration c = (7.2, 8.2). Dotted lines show density profiles for 
the corresponding fractions of most bound DM particles at z = 1 (assuming 
an isotropic NFW distribution function with virial radius and concentration 
given by the N-body halo of each galaxy), normalized to the same stellar 
mass. Solid lines show exponential profiles with the same amplitude and 
half mass radius as the dotted lines. 



2.3 Tagging criterion and the / mb parameter 

The particles we select for tagging are supposed to approximate 
the phase space distribution of the stars immediately after they 
form. Stars are the end result of dissipative collapse, so a basic re- 
quirement is that particles tagged with newly-formed stars should 
be deeply embedded in the potential well of their dark halo when 
we tag them. We achieve this by ranking DM particles in the halo 
by their binding energy and selecting all those more bound than a 
threshold value, corresponding to a fixed fraction of the mass of 
the halo. Following C10, we call this free parameter of the method 
the 'most-bound fraction', / m b- A value of / m b = 0.01 means we 
selected the 1 per cent most-bound particles. 

We stress that our model for the structure of merger remnants 
is not purely collisionless, because the Gil model explicitly in- 
cludes enhanced dissipative star formation (in the bulge compo- 
nent) during mergers. This is important because hydrodynamical 
simulations of galaxy mergers have shown that nuclear starbursts 
increase the central phase space density of merger remnants (Hern- 
quist, Spergel & Heyl 1993; Robertson et al. 2006; Hopkins et al. 
2008). We include stars formed in this way in our tagging, although 
we do not distinguish them from stars formed in the normal 'quies- 
cent' mode (in particular, we use the same / m b in starbursts). 

The choice of / m b is more-or-less arbitrary, but this freedom 
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allows us to tune the scale length of the in situ components of our 
galaxies in a predictable way. This is because, in an NFW potential, 
the surface density profile of dark matter more bound than a given 
energy is roughly exponential (at least for / m b < 10 per cent), with 
a scale radius that depends on the threshold energy. This result can 
be verified easily by integrating the cumulative energy distribution 
of an NFW halo up to a given fraction of its virial mass, and con- 
structing the corresponding density profile from the phase space 
distribution function. We have done this using numerical approxi- 
mations for the distribution function (and density of states) given by 
Widrow (2000) for a spherical NFW halo with an isotropic velocity 
distribution (more details are given in appendix |B|>. 

To illustrate this point, Fig.[2]shows the profile of in situ stars 
in two 'Milky Way' mass haloes from Millennium II (top and bot- 
tom panels), according to our full particle tagging model (dots) with 
/mb = 0.01 (blue) and / m b = 0.1 (red). Dotted lines show the 
profile we obtain using the Widrow (2000) distribution function 
to select the equivalent most-bound mass fraction at z = 1 when 
most of the stars in these galaxies form (the central regions of the 
are very stable, e.g. Wang et al. 2011). The dotted profiles are not 
exactly exponential because our procedure obviously imposes an 
energy threshold, which corresponds to a truncation radius. Solid 
lines show exponential profiles that have the same scale radius as 
the dotted profiles - these roughly approximate the diffusion of 
tagged particles across the initial energy threshold over time. Note 
that because we perform our tagging procedure at every snapshot, 
each new population in our full model will have a different amount 
of time to diffuse away from its initial configuration. 



2.4 Constraints on / m b from the galaxy mass-size relation 

From the above we conclude that our analytic approximation can 
reproduce the z — density profiles of particles tagged to rep- 
resent in situ stars with reasonable accuracy. This provides an in- 
tuitive understanding of why the parameter / m b sets the sizes of 
galaxies dominated by in situ star formation, and how those sizes 
vary with the scale of the dark matter potential. In this section, we 
use this understanding to determine a range of suitable / m b values 
empirically, by comparing the model M^—Rso relation to observa- 
tions of galaxies believed to be dominated by in situ star formation. 
By constraining / m b with reference only to in situ stars, the dis- 
tribution of accreted stars remains a valid prediction of our model. 
With a similar approach, C10 found that / m b = 0.01 gave reason- 
able agreement between their simulations and the Ah-R^o relation 
of dwarf satellite galaxies in the Local Group. We re-examine the 
choice of / m b because C10 considered only galaxies that were pre- 
dominantly satellites, with very different stellar and dark matter 
mass scales to those in our simulation. 

In Fig. [3] we show the median M+-R50 relation for all galax- 
ies in our full particle tagging model at z = (red, green and blue 
points, corresponding to / m b = 0.1,0.05 and 0.01 respectively). 
Each galaxy in the model contributes three values of R 5 q from three 
orthogonal projections to its M* bin. Model galaxies have a con- 
siderable scatter in R50 at fixed stellar mass (the 16-84 percentile 
range for / m b = 1 per cent is shown by the error bars on the blue 
points, other values of / m b have very similar scatter). 

At M+ < 10 11 , where galaxies in our model are dominated by 
in situ stars, we can use our analytic approximation to extrapolate 
our results relation below the limit of M* ^ 5 x 1O 1O M0 we im- 
posed when selecting galaxy merger trees for tagging. We take all 
central galaxies at z = in the Gl 1 Millennium II catalogue with 
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Figure 3. Points plot stellar mass Af* against projected half-mass radius 
i?50 for all galaxies in our full particle tagging models; colours correspond 
to different values of / m b ■ Cyan points show the masses and sizes of the 
z = galaxy sample at z = 1 for / m b = 1%- Dotted lines of the same 
colours (plus a magenta line for / m b = 3%) show the predictions of our 
analytic approximation for all central galaxies in the Millennium II simula- 
tion, using their in situ stellar mass and z = 1 halo properties. The dashed 
grey line and shaded region plot the late-type galaxy relation of Shen et al. 
(2003) and its lcr range, and the dash-dotted grey line the early-type re- 
lation of Guo et al. (2009). The solid cyan line shows the 2 = 1 relation 
for all galaxies from Ichikawa et al. (2012). Black points with error bars 
show the results of Kauffmann et al. (2003a) and their la range, also for 
the entire galaxy population. 



A'h ^ 1O S M0 and use the virial radius 0"2oo) and concentration^] 
c of their dark matter haloes at 2 = 1 to predict R 5 q at z = as de- 
scribed in appendix [B] These predictions are shown as dotted lines 
in Fig. [5] They agree roughly with the results of the full tagging 
model in the range where they overlap, although we find that R50 
is underestimated by ~ 33% for / m b = 1%. This disagreement 
is not surprising given the many approximations involved (includ- 
ing scatter in the age of the in situ population due to continual star 
formation and the crude representation of the central dark matter 
potential in our approximation). 

This extrapolation confirms that our tagging model produces 
a curved log 10 M*-log 10 -R50 relation similar to the observed re- 
lation of Shen et al. (2003, grey dashed line with la scatter). At 
lower A'h the model with / m b = 1% underpredicts the observed 
relation. Our approximation predicts that a model with / m b ~ 3% 
(magenta line) would be closer to the data; / m b = 5% is also plau- 
sible, as the Shen et al. (2003) relation may underpredict the sizes 
of edge-on galaxies by ~ 0.15 dex (e.g. Dutton et al. 2007). On 
the other hand, the Shen et al. (2003) relation we have plotted is 
that for 'late type' galaxies only (defined by Sersic index n < 2.5). 
Early-type galaxies with Ah < 10 11 Mq are significantly more 



4 c is not measured directly from the simulation; we use the result for an 
NFW profile that c = 2.16 »'20o/ r 'max where r max is the radius of maxi- 
mum circular velocity (e.g Cole & Lacey 1996). 
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Figure 4. Points with error bars show the median mass-size relation for el- 
liptical galaxies (Gil B/T > 0.9) predicted by our models with / m b=i% 
(blue) and / m b=5% (g reen )- Error bars indicate the 16-84 per cent range 
of .R50 in each mass bin. For log 10 M+/Mq > 11.5 we plot every galaxy 
in our / m b = 1% sample (blue triangles). We compare these points with 
individual galaxies from Guo et al. (2009) (grey points) and their best-fit 
linear relation (red dot-dash line; both corrected to a Chabrier IMF). Also 
shown are a stack of 14 nearby galaxies by van Dokkum et al. (2010) (light 
green triangle, a stack of 14 galaxies) and the relation of Hyde & Bernardi 
(2009) (green dashed line and 16-84 per cent range). 



compact (e.g. Shen et al. 2003). Because we are not able to model 
the morphology of in situ stars (notably disk structure) in a realis- 
tic way, it is difficult to reproduce observational classifications of 
early and late type galaxies, so it is preferable to compare to data 
that includes all galaxies. Kauffmann et al. (2003a) measured the 
distribution of half-light radius in the SDSS z band in bins of M* 
without making a morphological cut; their results (black squares) 
support 1% < f mh < 3%. 

Therefore, in order to bracket the range of fmb consistent with 
observations and illustrate the effect of this choice on our conclu- 
sions, we will use three fiducial values of / m b = 0.01, / m b = 0.05 
and fmh = 0.1 for the remainder of the paper. Within this bracket, 
we prefer not to fine-tune by fixing a 'best' choice of / m b, al- 
though values from 1% to 5% seem more likely. The simplicity of 
the model means that any choice must be approximate. Within our 
fiducial range, many other systematic effects beyond f m b (includ- 
ing the treatment of in situ star formation in Gl 1 and the cosmol- 
ogy of Millennium II) may determine how well our tagging model 
matches the observations (which themselves may not be accurate 
at the 0.1 dex level in R$o and M* corresponding to 1% changes 
in fmb, particularly in the conversion of flux to M*). 

We also note that in our model, the sizes of in situ stars evolve 
along the z — relation after z ~ 1, with very little change in 
the median i?so at fixed M* on average. This is illustrated by cyan 
points in Fig. [3] which show the z = 1 relation defined by the same 
galaxies that contribute to the blue points (i.e. our sample of z = 
central galaxies with fmb = 1%). This evolution agrees with the 
analysis of galaxy sizes in the GOODS North field by Ichikawa, 
Kajisawa & Akhlaghi (2012) who found very little change in R50 



at fixed M* between z — and z = 1 (their relation at z = 1 is 
shown by a cyan line). Early-type galaxies in the same mass range 
may undergo greater size evolution (e.g. Williams et al. 2010). 

Having determined a plausible range of / m b with reference 
to in situ stars, we now examine the predictions of the model for 
more massive galaxies that are dominated by accreted stars. In 
Fig. [3] the mass-size relation clearly steepens at M* > lO n M0. 
In this regime the / m b = 1% model follows approximately the re- 
lation for early-type galaxies found by (Guo et al. 2009) (which we 
have corrected from a Kroupa to a Chabrier IMF with a shift of 
Alog lo M*/M = -0.04), while f mh = 5% overpredicts R 50 
by ~ 0.15 dex and 10% by ~ 0.3 dex. Thus the model predic- 
tion i?5o in this mass range is slightly less sensitive to / m b than 
is for low mass, in-situ dominated galaxies (where R50 changes 
by ~ 0.5 dex over our fiducial range in / m b). In fact, as we will 
show in the following section, this change is still mostly driven by 
the strong effect of fmb on the scale of in situ stars even in very 
massive galaxies. The effects on the scale of the accreted compo- 
nent are even weaker, the most notable being an increase in scale 
at larger values fmb because more extended satellite galaxies are 
more easily stripped. 

Although a model with fmb ~ 5% gives an empirically rea- 
sonable scale for in situ stars in lower-mass galaxies, it appears to 
systematically overpredict the sizes of the most massive. To explore 
this apparent tension further, Fig.|4]plots the M*-7?5o relation for 
'elliptical' galaxies only (defined by Gil B/T > 0.9). The statis- 
tics of our simulated sample are poor above M* = lO n ' 5 M0: 
we plot each galaxy projection from the model as an individ- 
ual data point (blue triangle) in this region. These data are com- 
pared to observed relations for early-type galaxies from Hyde & 
Bernardi (2009) and Guo et al. (2009). Such comparisons are not 
straightforward, because of the variety of sample selections and 
measurement techniques for both in the literature. Furthermore, at 
M* > 10 11 ' 5 Mq where the disagreement of our models and the 
observational relations is most severe, there is an important system- 
atic uncertainty in the measurement of both M* and R50 in obser- 
vations which are not sufficiently deep to recover all the light in the 
outermost low surface brightness regions of very extended galax- 
ies. This can be a considerable fraction of the total light in galaxies 
with high Sersic index. Although in principle the bias due to this 
'missing' light can be corrected for (if the true surface brightness 
profile can be inferred), but it can also be compounded by other 
systematic effects in sample selection, background subtraction and 
profile fitting (Graham et al. 2005; Lauer et al. 2007; Blanton et al. 
201 1; Bernardi et al. 2012; Mosleh, Williams & Franx 2013). 

Hyde & Bernardi (2009) attempted to correct their relation 
(green dashed line and shaded la region) for the missing light bias 
within the limitations of the raw SDSS data. Their results agree 
with other SDSS-based relations, more so at lower Af* (e.g. Shen 
et al. 2003; Gadotti 2009; Williams et al. 2010). They also agree 
with the HST study of Auger et al. (2010) at higher M*. However, 
at high Mi, the Hyde & Bernardi (2009) relation does not agree 
with the somewhat steeper relation of Guo et al. (2009), who also 
analysed SDSS images. 

The key differences between Guo et al. (2009) and similar 
studies are that the authors preselect central galaxies, carry out 
their own background subtraction, and fit Sersic models to the sur- 
face brightness distribution, importantly with the Sersic index as a 
free parameter. This may allow them to more accurately recover the 
sizes of galaxies with very high Sersic index. We also plot the re- 
sults of van Dokkum et al. (2010) (green triangle), who stacked 
Sersic fits to individual deep images of 14 galaxies with mean 
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log 10 M*/Mq = 11.45 from the mass-selected and approximately 
volume-limited sample of nearby early types described by Tal et al. 
(2009, the OBEY survey). This point supports the Guo et al. (2009) 
relation. We also note that the datapoints plotted in figure 8 of Guo 
et al. (2009) for A'L > 10 11,5 M o (reproduced as grey points in 
Fig. [4]( lie systematically above the linear relation defined by the 
bulk of their early-type sample. This suggests that the relation may 
in fact curve upwards at higher M t (as claimed for example by 
Bernardi et al. 2011). 

The disagreement between different observations at high Ah 
suggests the uncertainty caused by undetected light in low sur- 
face brightness regions around massive galaxies may be substan- 
tial. With this in mind, we conclude that the median Rso of accre- 
tion dominated galaxies in our models are broadly consistent with 
observations provided 1% < / m b < 5%. In particular, the relation 
for simulated 'early type' galaxies at high A'h is steeper than that 
for in situ stars, in qualitative agreement with (Shen et al. 2003) 
and Hyde & Bernardi (2009). This is an improvement over the Gil 
semi-analytic model alone, which predicts an almost flat mass-size 
relation in this range of Af*. However, we also note that our rela- 
tion for early types does not agree with the observational consen- 
sus that ellipticals with M* < 10 are significantly more compact 
than late type galaxies with the same M* (Shen et al. 2003). On the 
other hand, because both late and early type galaxies are common 
below M* = 10 1 Mq, this result is sensitive to how early type 
galaxies are selected (and, perhaps, to the separation of central and 
satellite galaxies). 

2.5 Limitations of the method 

The most important limitation of our method is that we do not 
self-consistently model the gravitational effects of concentrating 
baryons in the cores of dark matter haloes. Our model is clearly 
inadequate wherever dark halo potentials are significantly modified 
by baryons. The work of C10 was restricted to the dwarf satellites 
of Milky Way-like galaxies, which (at least in the case of the MW) 
are highly dark matter dominated. Here, however, we tag dark mat- 
ter particles to represent stars in galaxies that are at least as mas- 
sive as the Milky Way itself. We cannot model the dynamics that 
govern the formation and secular evolution of discs, adiabatic con- 
traction or the destruction of cusps in dark haloes by feedback, all 
of which are thought to be important in the cores of massive galax- 
ies (Navarro, Eke & Frenk 1996; Gnedin et al. 2004; Dutton et al. 
2007). We proceed with this very simple model nevertheless as a 
first (if approximate) step towards modelling the spatial distribution 
and other properties of the diffuse stellar component of galaxies in 
the context of a realistic model of galaxy formation. Our model for 
in situ stars is only intended to serve as a means of creating initial 
conditions for accretion events with roughly the right scale and con- 
centration, and as a means of quantifying the relative contribution 
of accreted stars at different radii. 

Several other effects that could be important in real galaxies 
are not included in our model. If we were to include the gravity 
of stars in our central galaxies in a self consistent way, this would 
of course alter the overall structure of the potential and the regions 
of phase space accessible to accreted and in situ stars. Importantly, 
it could also disrupt satellites more easily, through tidal shocking 
and enhanced (stellar) dynamical friction. In our model, satellites 
passing near the centres of their hosts do not suffer these effects, 
which artificially favours their survival. Neglecting these processes 
has the opposite effect to the numerical limitations of our model 
(see below) which tend to enhance the rate of satellite disruption. 



Our tagging method cannot model rotationally supported discs 
because we tag particles based only on energy, and in any case 
there are very few coplanar dark matter particles with appropri- 
ate angular momenta in the simulation. Therefore even late type 
galaxies in the Gl 1 models are represented as dispersion-supported 
spheroidal systems in our particle tagging Although, as discussed 
in the previous section, their surface density profiles are exponen- 
tial and their scale lengths can be roughly matched to observations 
by choosing / m b appropriately, the energy and angular momentum 
distribution of stars stripped from these galaxies may be unrealis- 
tic. Finally, we also neglect the possibility that in situ stars form 
on very loosely bound orbits far away from galaxies, for example 
from cold gas clumps stripped from satellites or ejected in galactic 
fountains. 



3 RESULTS FOR INDIVIDUAL GALAXIES 

Before we present surface brightness profiles of galaxies averaged 
in bins of stellar mass and halo mass, we illustrate the basic out- 
put of our model with 12 examples of individual galaxies from our 
sample of 1872. 

3.1 Images 

Fig.|5]shows 2D images of projected stellar mass surface density for 
twelve of our simulated galaxies. In these images the stellar mass 
associated with each tagged dark matter particle has been smoothed 
over a conventional 'SPH' cubic spline kernel (e.g. Springel 2005) 
with a projected scale radius equal to the 3D distance to the 64th 
nearest dark matter neighbour (counting all particles, not just those 
with stellar populations). Satellite galaxies (tagged particles in self- 
bound subhalos) have been removed, so the remaining 'lumps' on 
small scales are due to the shot noise of the particle distribution. 

The galaxies labelled A-D in Fig. [5] represent the lowest 
halo masses in our sample, with A/200 ~ 10 12 Mq and M t ~ 
6 x 10 10 M (similar to the Milky Way, e.g. ; Watkins, Evans & 
An 2010; McMillan 2011). Those labelled E-H are more massive 
isolated haloes with M 20 o « 10 12,5 M and slightly larger M 4 , 
while those labelled I-L represent the central galaxies of groups 
and poor clusters, with 10 1325 < M 2 oo < 10 14 M Q and M* up to 
2 x 10 11 M . The Gil model predicts that B, C, D, G and K have 
B/T < 0.2, and J has B/T ~ 0.4. The other examples fall into 
our 'early type' category denned by B/T > 0.9. 

We have chosen a colour scale for surface mass density to 
illustrate the observability of different features. Regions with red 
colours are readily observable using conventional imaging tech- 
niques (for example in the Sloan Digital Sky Survey; [i r < 
25 mag arcsec -2 ) whereas those with yellow/green colours re- 
quire careful deep imaging (25 < fj,y < 28 mag arcsec" 2 , e.g. 
Martfnez-Delgado et al. 2010b). Careful reduction of SDSS im- 
ages (particular those in Stripe 82) can achieve 26 < /J, r < 
28 mag arcsec -2 (e.g. Kaviraj 2010; Bakos & Trujillo 2012). The 
next generation of imaging surveys will reach these depths rou- 
tinely; the Large Synoptic Survey Telescope (LSST) is expected to 
approach /iy ~ 29 mag arcsec -2 over a large fraction of the sky 
after 10 years (LSST Science Collaborations et al. 2009). 

Regions in Fig. [5] with blue colours are extremely challenging 

5 We can still separate galaxies by morphology using the predictions of the 
Gil model, which are independent of our tagging approach. 
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Figure 5. Projected stellar mass surface density in regions 500 X 500 kpc around 12 simulated galaxies. In each of the three columns, M200 and M* 
are approximately constant, representing (from left to right) Milky Way-like galaxies, massive isolated galaxies, and group/cluster central galaxies. These 
examples are referred to in the text by their labels A-L (blue where Gil B/T ^ 0.9 and red where Gil B/T > 0.9). More details are shown in Fig. [6] 
Satellite galaxies are included in the model but are not shown for clarity. £ = 7.0 (6.0, 5.0) Mq kpc -2 corresponds to a V band surface brightness of 
~ 24.8 (27.3, 29.8) mag arcsec -2 assuming a mass-to-light ratio of 2.5 (see Fig.|6]for this approximate conversion scale). We have imposed a minimum 
surface density of £ = 3.0 Mq kpc -2 (~ 35 mag arcsec -2 ) with Poisson noise to create the 'background' in these images. The lumpy appearance of the 
diffuse light is due to shot noise in the adaptively smoothed particle distribution. 
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Figure 6. Solid lines show stellar mass surface density profiles for the galaxies shown in Fig. [5] (identified by letters A-L). Line colours correspond to Gil 
B/T values as shown in the legend. Short black solid lines show slopes of -1.5, -2.0 and -2.5 for reference. The thick dashed black lines are the average total 
density profiles for galaxies with the same halo mass as the examples, discussed in section [4] Dotted lines show the density of in situ stars only. A vertical 
dashed line marks the effective force softening scale. 



to observe in individual galaxies. Nevertheless, for small areas and 
small numbers of nearby galaxies, surface densities equivalent to 
/j, r ~ 31 mag arcsec -2 have already been explored by resolved 
star surveys (e.g. Bland-Hawthorn et al. 2005; Barker et al. 2009; 
Bailin et al. 201 1 ; Radburn-Smith et al. 201 1). Supernovae, plane- 
tary nebulae and globular clusters are expected to trace this under- 
lying stellar distribution. Perhaps the most promising method for 
direct studies of these faint regions is stacking the surface bright- 
ness profiles of many similar galaxies (Zibetti et al. 2004, 2005; Tal 

6 vanDokkum 2011). 

3.2 Individual density profiles 

The solid lines in Fig. [6] show the stellar mass surface density pro- 
files of the 12 galaxies from Fig. [5] Galaxies in the three panels 
(representing different ranges of M200 and M t ) show clear differ- 
ences in the shape of their profiles, as well as the expected differ- 
ences in amplitude between galaxies of different mass. Whereas 
there are no clear differences in the appearance of the galaxies 
in the first two columns of Fig. [5] the late type galaxies (B, C, 
D and G) clearly show two structural components in Fig. [6] an 
exponentially-declining inner profile that breaks to a shallower 
slope at radii between 10 kpc and 30 kpc. This break is not seen 
in the profiles of their early type counterparts (A, E, F, H), which 
have roughly constant slope from 5 kpc to 100 kpcj^] 

The inner regions of these galaxies (< 10 kpc) are very similar 
to each other, even for different halo masses and at radii well out- 
side the effective force softening scale (shown by a vertical dashed 
. Coloured dashed lines in Fig. [6] show that in situ stars dom- 
inate these regions. Although the in situ component grows in size 

e We stress that these profiles and the classification of the galaxy as early 
or late type in the Gil model are not directly connected; they are only 
correlated by their individual dependence on the mass accretion history. 

7 Forces are exactly Newtonian for particle separations greater than 2.8e = 



with halo mass, the mass of the accreted component that dominates 
at larger radii increases even more rapidly. The surface density at 
~ 100 kpc, which is almost entirely contributed by accreted stars 
in all our examples, increases by three orders of magnitude across 
the panels. 

Therefore, the physical origin of the break mentioned above 
is the transition from regions dominated by in situ stars to regions 
dominated by accreted stars. This implies a clear connection be- 
tween the typical shapes of galaxy surface brightness profiles and 
the relative fraction and distribution of accreted stars. The rest of 
this paper will focus on the importance of accreted stars in chang- 
ing the surface brightness profile shape as a function of halo mass, 
using all the galaxies in our sample. 

3.3 Stellar mass-to-light ratios 

For simplicity we present most of our results in terms of stellar 
mass surface density, and only convert to surface brightness where 
necessary using fiducial, constant values for stellar mass-to-light 
ratio in a given band, M+/L. This is motivated by Fig. [7] where 
were show the results of applying stellar population synthesis to 
each of our tagged particles, based on the age and metallicity of 
their stars. For the sake of illustration, we do this for three exam- 
ple galaxies chosen to be typical of galaxies of similar mass (it is 
possible to find individual galaxies with quite different stellar pop- 
ulation gradients in the model). Following Gl 1, we use the Bruzual 
& Chariot (2003) SED templates and a Chabrier (2003) IMlQ 

3.84 kpc with h = 0.73 and Plummer-equivalent softening length e = 
1 hT 1 kpc fixed in comoving coordinates. 

8 Because our stellar tagging is done in postprocessing, the ages of each 
population are discretized by the output times of Millennium II. When cal- 
culating ages, we assume stars are formed half-way between snapshots. 
These snapshots are spaced ~ 300 Myr apart from z ~ 1 to z = 0, 
which is not sufficient to model the rapid changes in flux (mostly at short 
wavelengths) from very recent star formation. However, very young pop- 
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Figure 7. Population synthesis results for three galaxies with (from left to right) log 10 M20o/MQ = (12.0,12.4,14.2) and log lo M*/M0 = 
(10.8, 11.0, 11.3). Bottom row: surface brightness profiles in SDSS ugri bands and the Johnson V band (solid lines) and approximations assuming con- 
stant stellar mass-to-light ratios given in the text (dots). A vertical grey dashed line marks the effective force softening scale. Middle row: colour gradients; for 
reference dotted lines (the same in each panel) show fiducial gradients of A u _ r = —0.2 (e.g. Peletier et al. 1990; Tal & van Dokkum 201 1), A 3 _; = —0.09 
and A 9 _ r = —0.07 La Barbera et al. 2011). Filled circles show the mean colours predicted by the Gil model (without dust correction). Top row: pop- 
ulation synthesis stellar mass-to-light ratio as a function of radius corresponding to the surface brightness profiles in the lower panels; dotted lines indicate 
Nh/L = 2,3,4,5. 



We draw two conclusions from this figure. First, over the 
mass range of our sample, converting stellar mass surface den- 
sity to surface brightness assuming a constant M*/L for a given 
band is a reasonable approximation to the results of our full stellar 
population modelling, especially in the outer parts of our galax- 
ies (within a factor of 2 from 10 to 100 kpc). The dots overplot- 
ting the profiles in Fig. [7] (solid lines) correspond to A/*/ ' L ugr i — 
18.0, 4.0, 2.0, 1.5 and M„/L v = 2.5. (In the same way we find 
approximate values of M*/Lb,r,i = 6.0,2.0,1.5; for clarity 
these lines are not shown in Fig.^J. Only the inner core of the low- 
est mass galaxy in these three examples has a significantly different 
M^/ L from the outskirts. 

Second, consistent with the behaviour of M*/L, the colour 
profiles of our galaxies are generally very flat. The principal reason 
for this is that in the Gil model, recent star formation in haloes 
with A/200 > 10 13 Mq is strongly suppressed by AGN feedback 
(Kauffmann & Haehnelt 2000; Croton et al. 2006; Bower et al. 
2006). The central galaxies of such haloes typically grow through 
mergers with 1-2 objects of similar stellar mass after z ~ 1. These 



ulations only make a small contribution to the masses of galaxies in our 
sample. The underlying semi-analytic model of Gl 1 has better time resolu- 
tion, because it calculates star formation on a grid of 20 'sub-steps' between 
the simulation snapshots. 



have old and metal rich stellar populations very similar to those of 
the central galaxy because their star formation has also been sup- 
pressed by AGN feedback (see Section|7]l. Also, as seen in Fig. [6] 
and discussed further in Section [4] the density profiles of the ac- 
creted and in situ components are more similar for more massive 
galaxies, which also reduces stellar population gradients. Colour 
gradients are stronger at lower M200 and in bluer colours, as in 
our lowest mass example, where the accreted stellar halo is con- 
siderably older, more metal poor, more extended and less massive 
than the recently-formed in situ component. Our models are gen- 
erally consistent with the observations of Peletier et al. (1990), La 
Barbera et al. (201 1) and Gonzalez-Perez, Castander & Kauffmann 
(2011), which we overplot as dashed lines in the middle row of 
Fig. [7] and with the mean colours predicted by the Gil model for 
the same galaxies (without dust correction; dots at 10 kpc in the 
middle panel). 



4 AVERAGE SURFACE DENSITY PROFILES 

Figure[8]shows the main results of this paper: the median profiles of 
stellar mass surface density in circular annuli, for all the galaxies in 
our sample in logarithmic bins of 0.5 dex in dark halo virial mass 
(Af2oo)- In this figure we only show results for / m b = 0.01, for 
clarity. We do not take projections along the principal axes of the 




Figure 8. Median profiles of circularly averaged stellar mass surface density, for accreted stars (red dashed lines) and in situ stars (red dotted lines), in five 
logarithmic bins of dark halo virial mass (11.5 < log 10 M200 /Mq < 14.0, shown in the top right of each panel along with the number of galaxies in each 
bin). A blue solid line shows the median profile (for / m b = 0.01) combining the accreted and in situ components; a light blue region around this line indicates 
the 10-90 per cent scatter of the median profile. (Every halo contributes three orthogonal projections to each median.) Arrows indicate the half-mass radius of 
the median profiles (from left to right, in situ stars, all stars and accreted stars). Grey lines (dotted, dashed and solid) reproduce the corresponding red and blue 
lines from the 12.5 < log 10 M200/MQ < 13.0 panel. A purple line and pink shading show the median dark matter density profile and the corresponding 
10th to 90th percentile region. These DM profiles are not measured directly from the simulation but approximated by assuming projected NFW profiles with 
concentration c = r v ; r /r s derived from the ratio of the virial radius and radius of maximum circular velocity using the relation c = 2.16 -R200 / Kmai (e.g. 
Cole & Lacey 1996). A black horizontal bar shows the range of the virial radius i?200 m eacn mass bin, and a vertical dotted black line indicates the effective 
softening scale 2.8e. The scale on the right of the lower central panel gives an approximate conversion from to surface brightness (in Vega magnitudes per 
square arcsecond) for the Johnson-Cousins V band, assuming Ty = M+/Ly = 2.5 [for alternative mass-to-light ratios, T' v , add 2.5 log 10 (T^/2. 5) to 
these values]. 



galaxies or align them in any other way, so the relative orientations 
of the profiles we combine are random. 

The profiles of total stellar mass surface density (blue lines) 
show a tight correlation in shape and amplitude with M200, with 
80 per cent of profiles differing by no more than 0.5 dex from the 
median at all radii in all halo mass bins (a large part of this scatter, 
shown by the light blue band, is due to the 0.5 dex width of our halo 
mass bins). The separate contributions of accreted and in situ stars 
explain the shape of the profiles and their variation with halo mass. 
In situ stars (red dotted lines) dominate the high surface brightness 
regions in haloes up to ~ 1O 13 M0, while accreted stars (red dashed 
lines) dominate at radii greater than 30 kpc and in more massive 
haloes. Regardless of M200, neither the two subcomponents nor 
the combined stellar profile follow the NFW distribution of the dark 
matter (purple lines). 

In Section |3.2| we discussed an interesting relationship be- 
tween the bulge-to-total stellar mass ratio B/T (predicted by the 
Gil model) and the shape of the surface density profiles we obtain 



with particle tagging. In Figure [9] we show that this result holds 
on average, plotting the / m b = 1% profile in the bin 12.0 < 
log 10 M200 < 12.5 (top-centre panel of Figure [8J separately for 
galaxies with B/T > 0.9 and B/T < 0.9. The latter profile shows 
a clear inflection at ~ 10 kpc, which corresponds to the radius at 
which accreted stars begin to dominate over in situ stars. The ab- 
sence of this break in galaxies with B/T > 0.9 reflects a greater 
contribution of accreted stars to the inner parts of galaxies and the 
action of violent relaxation in major mergers (the primary cause of 
high B/T in the Gl 1 model), which makes the profiles of the two 
components more similar. Figure [9] shows that profiles of galaxies 
with low and high Gil B/T can still be distinguished if we use 
a much higher value of / m b, although the accreted-in situ break is 
then much less clear. We caution that the value of B/T is calcu- 
lated in a simple manner and is influenced by other processes in the 
Gil model (such as disk instabilities and starbursts). In individ- 
ual cases it does not correspond perfectly to a separation between 
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Figure 9. Average stellar mass surface density profiles for / mb = 1% 
(solid lines) and 10% (dot-dashed lines) in the halo mass bin 12.0 < 
log 10 M 2 oo < 12.5, subdivided into galaxies with B/T < 0.9 (N = 
704, blue lines) and with B/T > 0.9 (N = 132, red lines). 



two-component and quasi-power-law profiles in our tagged particle 
model. 

We find that the accreted and in situ surface density profiles are 
reasonably well described by Sersic functions in all haloes. How- 
ever, a single Sersic model is only appropriate for the total profile 
in halos more massive than M 2 oo ~ 10 13 Mq, where the accreted 
component dominates at all radii. Even in these haloes, the in situ 
component makes a significant contribution, and a slight change 
in shape due to the transition between accreted and in situ stars is 
apparent at R < 10 kpc. Thus a double-Sersic profile provides a 
good description of the total stellar mass surface density profile in 
all haloes. 

Fig. [TO] shows how the effective (i.e. projected half-mass) ra- 
dius R50, amplitude £50 and index, n, of Sersic model fits to the 
surface density of accreted stars (left) and in situ stars (right) vary 
with A/200- We further subdivide galaxies by their B/T ratio in 
the Gil model. Table [T] presents linear fits to these trend^] The 
scaling behaviour of the in situ stars is largely a consequence of 
our particle tagging method (as described in Section \2.3\ , so we 
focus on the scaling of the accreted component. These scalings are 
not easy to predict without direct simulation because they depend 
on the masses and star formation histories of all the progenitors of 
a given halo, as well as the collisionless dynamics of dark matter 
structure formation. 

Looking first at 'elliptical' galaxies with B/T > 0.9, the 
effective radius of accreted stars scales (very roughly) as Rao oc 



M. 



2/3 



Changes to / mb scale profiles at all M200 by a similar 



amount, so this relation is not sensitive to the choice of / mb . The 
dependence of Rbo on M200 and / mb is reflected in the behaviour 



9 Some of the relations are clearly not linear; thus the fits are only a guide 
to the average behaviour and are less reliable as a way of reproducing our 
results than simply interpolating the data points we show in Fig. 1101 



of £50 shown in the middle panel (total stellar mass is roughly con- 
served with changes in / mb in the range we adopt). However, the 
slope of the relation between £50 and M200 also depends on / mb , 
with a considerably steeper relation at / mb = 0.01. This is because 
profiles in the models with / mb = 0.01 show an increase in n 
(i.e. higher concentrations) at lower Af 2 oo- This trend vanishes for 
fmb = 0. 1, for which n ~ 3 at all M^oo- The reason for this is not 
clear and might be better understood with controlled experiments 
like those of Hilz et al. (2012a). Furthermore, although the differ- 



10 



ence in Sersic index at M200 ~ 10 Mq seems large in Fig 
the difference in the actual profiles is smalf^] 

We conclude that (over the M 2 oo range of our study) ellipti- 
cal galaxies in more massive haloes have more extended accreted 
stellar mass surface density profiles, but the shape of these profiles 
is n ~ 3 in all haloes. Lowering the typical binding energy of pro- 
genitor galaxies results in accreted profiles that are more extended 
and slightly more concentrated. 

We now examine the 'late type' galaxies with B/T ^ 0.9. 
The effects of / mb on the accreted profile appear to be less in 
these galaxies, especially at low M 2 oo- Roughly, they have R50 oc 
M 2 5q 6 , a shallower relation than for the 'ellipticals'. However, lin- 
ear relations are not good fits to these trends and the highest- M200 
points (above M200 ~ 10 12 ' 5 Mq) show hints of steepening to a 
slope similar to that of the accreted stars. 

Regardless of this, it is clear that at all M200. the accreted 
haloes of late-type galaxies are much more extended than those 
of ellipticals. With marginal significance, they also have a very 
slightly lower n at fixed Mq,oo- The larger scale radius of the haloes 
of later type galaxies is likely to be due to the fact that they are built 
by several minor mergers which deposit stars on weakly-bound or- 
bits. In the case of ellipticals the distribution of accreted stars is 
more likely to reflect violent relaxation following a low mass ra- 
tio merger (which deposits a larger fraction of stars in the centre 
of the remnant). On this basis, we speculate that there may be a 
readily observable systematic difference in the scale radius of the 
outer envelopes of SO galaxies compared to ellipticals of similar 
total luminosity. 

The fact that n ~ 3 and R 5 q ~ 10 kpc are approximately con- 
stant for accreted stars in low-mass late-type galaxies suggests that 
accretion cannot produce the compact exponential bulges observed 
in many disk galaxies (at least in the M200 range we study). This 
is the case even though we include all the mergers in the history of 
each galaxy (including those between compact progenitors at high 
redshift). This supports the idea that compact exponential bulges in 
late type galaxies are formed in situ (i.e. through starbursts, perhaps 
in the same way as the cores of elliptical galaxies) or through the 
secular evolution of disks. 

Our linear approximations for the Sersic parameters of the 
two components are only reliable within the M200 range we study, 
and should not be compared with observations outside this range. 
At high M200, this is not an intrinsic limit of the particle tagging 
method, but rather because we only have a small number of haloes 
this massive in our volume. On the other hand, below the lower 
M200 limit of our sample, our results are limited by the resolu- 
tion of the simulation and our crude model for the in situ stars that 



10 In this regard, two points are worth noting; first, the fit to the trend in n 
is influenced by the lowest point (which is close to the resolution limit in 
-R50), and second, comparison to the constant-n case of / mb = 0.1 shows 
that most of the order-of-magnitude difference in E50 between / mb = 0.1 
and / mb = 0.1 can be explained by the change in R50. 
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Figure 10. Parameters of Sersic model fits to the median surface density profiles of accreted stars (left column) and in situ stars (right column) for simulated 
galaxies in narrow bins of M200 (containing 100 galaxies each). From top to bottom, rows show the three Sersic parameters as a function of M200: projected 
half mass radius R50, amplitude S50 and index n, respectively. The lines in each panel are linear fits to the data points, with colours corresponding to models 
with / m b = 0.01 (blue), 0.05 (green) and 0.1 (red). Each panel has two sets of lines and points, corresponding to galaxies with B/T ^ 0.9 (circles/solid 
lines) and B/T > 0.9 (crosses/dashed lines). Note that the first two rows/fits are log-log and the last is log-linear. A horizontal dotted line in the first row 
indicates the effective force softening scale. Error bars on the blue points in the top left panel are derived from a bootstrap resampling of the galaxies in each 
bin; the size of these error bars is similar for all other points, where they are not shown. 
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14.69 
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0.1607 
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-0.3434 
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0.2561 
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0.10 Acc 
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0.6262 
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B/T <: 0.9 
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-2.487 
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4.138 


0.05 Ins 
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Table 1. Linear approximations to the scaling of Sersic parameters with halo mass M200 derived from fits to our median surface density models. These are 
presented separately for ellipticals (Gl 1 B/T > 0.9, top) and late type galaxies (Gl 1 B/T < 0.9, bottom). From left to right columns show the most bound 
particle fraction / m b> the component (either accreted or in situ), and the slope a and intercept b for relations of the form y = C1M2QQ + ^> where y represents, 
respectively, log of the projected half-mass radius, log of the surface density at R50 and Sersic index. 




Figure 11. Surface density profiles in bins of total stellar mass, M*. Note that the final two bins are wider than the first three. Different line colours correspond 
to different choices of / m t, , which mainly affect the in situ component; the effect on the accreted component is negligible except in the most massive galaxies 
(see Appendix A). The light blue region indicates the 10-90 per cent scatter of the median profile for the / m t> = 0.01 case. Arrows show half mass radii and 
black lines show the range of -R200 m eacn bin- 
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dominate in low-mass haloes. We note that the simple semi-analytic 
approximation to our tagging method described in Section [2~3| pre- 
dicts that, in a higher resolution simulation, the value of n for in 
situ stars would asymptote to 0.5 < n < 1.0, log 10 E50 would de- 
crease as the result of the steep relation between M* and M200, and 
R50 would continue to decrease with Ahoo- This behaviour is in 
qualitative agreement with observations and, in turn, suggests that 
our conclusions regarding the accreted component (which is gen- 
erated by stripping lower mass haloes) are robust. This can only be 
properly tested with a higher resolution simulation, however (see 
Appendix [A i. 

Fig.^ repeats Fig. [5] but bins galaxies by stellar mass, Af*, 
rather than A/200 • The strong trends with A/200 are less clear in the 
case of A/*, except in the most massive stellar mass bin. The cause 
of this is the considerable scatter in the M+-M200 relation (e.g. 
Guo et al. 2010). Below their half-mass radii, the median surface 
brightness profiles are hard to distinguish across a range of 10.5 < 
log 10 M* < 11.5. A stronger evolution of the profiles with M* can 
be seen at larger radii, ~ 30-100 kpc. Fig. [TTJ also illustrates the 
effects of varying / m b, over the range 1-10 per cent. Relative to the 
1 per cent profiles (blue lines), the 10 per cent profiles (red lines) 
are ~ 30 per cent more extended and have lower central surface 
density (see section |2~3| >. This rescaling of the in situ component 
is the main effect of changing f m h. The most significant impact of 
higher / m b on the accreted component is in the very outer regions 
of the most massive galaxies, where slightly more stars are found 
on weakly bound orbits (perhaps as the result of earlier stripping). 
Again we conclude that our density profile results do not depend 
strongly on the value of / m b, at least within a range consistent with 
the observed size-mass relation. 



5 COMPARISON WITH LITERATURE DATA 

In Fig.[l2]we compare the surface density profiles shown in Fig.[8] 
(binned by A/200) with deep observational data from a variety of 
sources (summarised in Table [2}. We do this to illustrate the va- 
riety of different surface density profiles that have been reported 
in the literature, rather than to match any particular observation. 
In most cases our choice of an M200 bin for each observational 
dataset is based on the observed stellar mass (e.g. Guo et al. 2010; 
Moster et al. 2010) and is thus very rough. Where authors have 
presented their data in terms of stellar mass surface density, we use 
their values directly. Otherwise, we assume a galaxy-wide mass- 
to-light ratio appropriate to the observed photometry (as given in 
section [3^2] and listed in Table [2}. There are many systematic dif- 
ferences between these datasets (including photometric bandpass, 
surface brightness dimming corrections, A"-corrections, cosmology 
and in some cases, the choice of IMF). We have attempted to cor- 
rect for these differences where necessary. Such corrections amount 
to less than 0.1 dex in most cases. 

In the halo mass range 12.0 < log 10 A/200 /Mq < 12.5 
we show data from galaxies comparable to the Milky Way and 
M31 (the dark matter haloes of these galaxies may have A/200 ~ 
10 12 Mq, e.g. Watkins et al. 2010). The composite / band model 
of M31 from Courteau et al. (2011, cyan line) includes a power- 
law stellar halo of index 7 ~ —2.5, which implies a higher surface 
density at R > 100 kpc than our models predict. The Courteau 
et al. (201 1) halo is based on the data of Gilbert et al. (2009), who 
measured the surface density from star counts in individual fields 
(orange dots), some of which may contain substructure. At R < 
50 kpc these data agree well with the average profile of accreted 



stars in the model (dashed grey line). The galaxies M81 (Barker 
et al. 2009), NGC 2403 (M* ~ 1O 1O M ; Barker et al. 2012), NGC 
1087 (M^ ~ 10 10 ' 4 A/q) and NGC 7716 (A/* ~ 10 10 ' 5 A/ Q ; Bakos 
& Trujillo 2012) show similar profiles. All show hints of breaking 
to a shallower slope beyond 15-20 kpc, although these upturns oc- 
cur close to the limiting depth of the observations. 

In the same panel we compare with the GIMIC SPH simula- 
tions (Crain et al. 2009), the only other large cosmological simu- 
lation of stellar haloes in this M200 range (though see also Croft 
et al. 2009). Font et al. (201 1) stacked ~ 400 galaxies from GIMIC 
(with a particle mass 7.7x larger than Millennium II and the same 
force softening length). They find circularly averaged density pro- 
files that are well described by a concentrated in situ component 
and a diffuse accreted component, although the transition between 
the two components is less obvious in their profiles than in our 
/mb = 1% model. The shallower outer slope seen in the GIMIC 
simulations is in better agreement with the M31 data of Gilbert 
et al. (2009) at ~ 100 kpc than our models. 

In our simulations, accreted stars begin to dominate the sur- 
face brightness profile in the halo mass range 10 12 ' 5 < A/200 < 
1O 13 M0. In this range we compare with two galaxies from Bakos 



(M* 



1O 1O,9 M0). These profiles have higher amplitude at 



30 kpc relative to the galaxies in the previous mass bin. The sim- 
ulated galaxies show a similar increase, due to their much greater 
mass of accreted stars (dashed line). However, some features in the 
NGC 7761 and NGC 1068 profiles are not reproduced in the mod- 
els, specifically inflections at ~ 10 kpc and ~ 30 kpc and the rapid 
fall-off at R > 30 kpc. We note that Bakos & Trujillo (2012) report 
signs of tidal disturbance and substructure in these galaxies. 

In the range 13.0 < log 10 A/ /Mq < 13.5 we compare with 
the stacked surface brightness profile of ~ 42000 SDSS LRGs (lu- 
minous red galaxies, thought to be mostly group/cluster centrals; 
Mi, ~ 10 11 M©, A/200 ~ 10 13 ' 2 M Q ) at (z) ~ 0.34 from Tal & 
van Dokkum (2011). These data are a good match to the surface 
density profiles of our simulated galaxies from 10-100 kpc. Below 
10 kpc the simulated profile is steeper than the data, although this 
region is sensitive to our treatment of the in situ component. The 
simulation does not reproduce the upturn in the observed profile at 
~ 100 kpc, which Tal & van Dokkum (2011) attribute to limita- 
tions in their correction for residual light from unresolved compan- 
ion galaxies. Our models supports this interpretation because they 
do not predict a separate physical component with a shallow den- 
sity profile that could explain the upturn. We also show the stacked 
profiles of 14 nearby early type galaxies from the volume-limited 
OBEY survey (Tal et al. 2009). This stack includes all ellipticals in 
the survey with stellar mass log 10 A/^/Mq = 10.45 ± 0.15 (ap- 
pendix D of van Dokkum et al. 2010). These data match the Tal & 
van Dokkum (201 1) stack at small radii and do not show an excess 
over our simulated profile at R > 100 kpc, possibly because van 
Dokkum et al. (2010) stack 2D fits to the OBEY galaxies rather 
than stacking the images directly. 

Finally, the most massive haloes in our simulation have 
13.5 < log 10 A/2OO/M0 < 14.0. We compare these with indi- 
vidual BCG profiles from Kormendy et al. (2009, NGC 4486/M87) 
and Seigar et al. (2007, Abell cluster cD galaxies NGC 3551, NGC 
4874, NGC 6173, UGC 9799 and GIN 478). The data of Seigar 
et al. (2007) are represented by their 'double Sersic' fits. Simu- 
lated profiles in this halo mass range match these profiles well 
from 10-100 kpc, although some of the observed galaxies could 
belong to haloes more massive than the median of our bin (e.g. 
log 10 A/ 2O o/M = 15.1 for NGC 4874/Coma and 14.3 for UGC 
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Figure 12. As Fig. [8] The surface density profiles in each bin of halo mass are shown in grey, and compared with observational data (and other simulations) 
given in Table[2] The assignment of galaxies to halo mass bins is approximate. 

Table 2. Surface mass density profile data shown in Fig. |12| From left to right, columns give: the range of halo mass to which the data (or simulations) are 
compared (log 10 M200 /Mq, corresponding to panels in Fig. |12| ; the target galaxy or galaxies; the source of the data; the symbol or line style used in Fig. |12| 
the photometric bandpass of the data; the stellar mass-to-light ratio we have assumed (where the original authors do no present their results in terms of stellar 
mass surface density); and comments on the data (LTG: late type galaxy). 
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Figure 13. Symbols show the average stellar mass surface density profiles obtained from stacks of SDSS DR9 r band images as described in the text and 
appendixjc] assuming a constant stellar mass to light ratio (open circles) and, where significantly different, a colour-dependent M/L (open triangles). Error 
bars approximate 'la' of the distribution of uncertainty in the average profiles combining Poisson errors in flux measurement with the sample variance of the 
stack (Njjjig given in each panel shows the number of galaxies in the bin). Coloured lines (blue, green and red) show stacks made from our simulations as 
in Fig. |ll| but here binning galaxies by their Petrosian mass M pe t (see text). The lower central panel shows the four SDSS profiles only (colours indicate the 
central mass of each bin). 



9799/Abell 2052; Reiprich & Bohringer 2002). We also show the 
results of Zibetti et al. (2005) who stacked SDSS images for BCGs 
at (z) ~ 0.34. This is deeper than any of the individual profiles and 
shows a clear excess over our models beyond 100 kpc. We note that 
mass-richness relations suggest that the typical halo mass of their 
sample is M 20 o > 10 14 2 M (e.g. Rozo et al. 2009). 

Subject to the crude way in which we have assigned galax- 
ies to halo mass bins and our simplification of a constant mass-to- 
light ratio, the range of shapes, amplitudes and scales seen in our 
simulated surface density profiles for accreted stars appear to agree 
qualitatively with observations spanning a similar range of M* . The 
most serious potential disagreement is at radii > 100 kpc, where 
the models appear to underpredict the surface density observed in 
M31 and the clusters stacked by Zibetti et al.. 



6 COMPARISON WITH STACKED SDSS DATA 

It is more appropriate to compare the average galaxy surface den- 
sity profiles from our models with similar averages constructed 
from large galaxy samples (such as the stacking analyses of Tal 
& van Dokkum 2011 and Zibetti et al. 2005) than with individ- 
ual galaxies or small surveys as we did in the previous section. We 



have carried out our own simple stacking analysis of massive galax- 
ies observed by the Sloan Digital Sky Survey (SDSS) Data Release 
9 (DR9; Ahn et al. 2012) in bins of stellar mass (as given by the 
MPA-JHU Value-Added Cataloguj 1 " 1 "}. Our method for construct- 
ing stacked images is described in appendix [C] 

The resulting density profiles are shown as open circles in 
Fig. [13] split into four bins of stellar mass, each of which is ob- 
tained from a stack of Ndr9 galaxies as indicated. The panel la- 
belled "SDSS data only" summarises these four profiles, showing 
a clear shift in amplitude from the least to the most massive bin, 
out to the largest measured radius. Each panel assumes a constant r 
band mass-to-light ratio (the average of the MPA-JHU M/L r val- 
ues in the corresponding MPA-JHU mass bin, ranging from 2.3 
to 2.8 from the first to last bin; see appendix [Cj but this result 
holds even if the same M / L r is used for all panels, or if a colour- 
dependent M/L r relation is used ( )]Bell03. However, these results 
do not show any significant change in profile shape of the kind seen 
in previous figures (including observations of individual galaxies). 

Each panel compares our SDSS stacks to the average profiles 
of simulated galaxies binned by their Petrosian stellar mass M pGt 



11 http : / /www . mpa-garching . mpg . de/SDSS/DR7 
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Figure 14. Left: points show the total mass of accreted stars, -M acc as a function of Af200. split by bulge-to-total ratio as shown in the legend. Right: points 
show M acc as a function of stellar mass, M * . Grey lines correspond to fractions of the total stellar mass as indicated. 



(blue, green and red lines for / mc , = 0.01, 0.05 and 0.1, respec- 
tively). We use M pct here to reproduce the bias between the 'true' 
stellar mass of a galaxy and measurements based on 'corrected' 
aperture photometry with a finite surface brightness limit, such as 
the modelMag magnitudes from which the MPA-JHU masses are 
derived (see appendix [C}. The galaxies in each bin of Fig. |13|are 
therefore different from those in the corresponding bin in Fig. 1 1 1 1 
The observed and simulated profiles in Fig.[T3]should only be com- 
pared at R > 5 kpc, because the point spread function (PSF), 
which we have not deconvolved, dominates the observed profiles 
at smaller radii, and numerical softening affects the simulated pro- 
files. 

The lack of a clear trend in profile shape with M pct is also 
apparent in the simulations. This can be explained by two effects 
which mix massive galaxies with high Sersic index into lower stel- 
lar mass bins: scatter in the M*-M2oo relation and the systematic 
underestimation of M* by Af pct , which is more severe for more 
massive (high n) galaxies. Taking these effects into account, the 
data agree remarkably well with the simulations for / m b ~ 0.05. 
A slightly lower value of / mc , (> 0.02) would give even better 
agreement. 



7 ORIGIN AND STRUCTURE OF STELLAR HALOES 

Having demonstrated that our particle-tagging model produces sur- 
face density distributions that agree reasonably well with observa- 
tions, we now use it to examine the origin of stellar haloes as well 
as the relationships between the properties of central galaxy and the 
structure of their diffuse light. 



7.1 Mass of accreted stars 

The left panel of Fig. [14] shows M acc , the total mass of accreted 
stars in each of our galaxies, as a function of virial mass, M2oo- 
The right panel shows A/ acc as a function of the total stellar mass 
of the system, M*. A/ acc increases much less steeply with M200 
above a characteristic mass A/200 ~ 10 12 ' 5 M©. This corresponds 
roughly to the transition mass predicted by the Gil model, above 
and below which the galaxy mass function is dominated by early 
and late type galaxies respectively. We see in Fig.[l4]that elliptical 
(Gil B/T > 0.9) and late type galaxies have a clean separation 
at an accreted stellar mass of Af aC c ~ 2 x 1O 1O M , which (as 
shown in the right panel of Fig. 1 14} corresponds to ~ 30 per cent 
of the total (accreted plus in situ) central galaxy stellar mass M* 
(for all M200). This simply reflects the fact that the 'major merger' 
criterion for the destruction of disks (hence formation of elliptical 
galaxies) in Gl 1 is a progenitor mass ratio of 1 :3 or lower. 

Purcell et al. (2007) made similar predictions for Af acc as a 
function of A/200 at z — using prescriptions for halo assembly 
histories and subhalo orbital properties based on numerical simula- 
tions. As their results demonstrated, the relationship between A/ acc 
and A/200 seen in Fig.[l4]is a natural outcome of the CDM progen- 
itor halo mass function and the relation between A/* -A/200, which 
is thought to be roughly monotonic (e.g. White & Frenk 1991 ; Ben- 
son et al. 2000; van den Bosch, Yang & Mo 2003; Eke et al. 2004). 
To first order, this means that M acc is set by the typical ratio of 
A/200 between the few most massive accreted progenitors and the 
main halo (~ 0.1-1%). However, A/*-A/2oo relations derived from 
galaxy abundance matching (like those predicted by semi-analytic 
models) show an inflection corresponding to a peak in galaxy for- 
mation efficiency at A/ pca k ~ A/200 ~ 10 12 Mq (Guo et al. 2010; 
Moster et al. 2010, e.g.)). This creates two regimes in the scaling 
of A/ a cc with A/ 2 oo- 

In haloes with A/200 J; A/ pca k, in situ star formation effi- 
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Figure 15. Solid blue and red lines show the median ratio of accreted stel- 
lar mass (both halo and bulge) to the total galaxy mass (sum of accreted 
and in situ stars). Colours separate galaxies by B/T as given in the legend. 
Dashed lines of the same colour enclose 90 per cent of the distributions. 
Black solid and dashed lines correspond to the same quantities for the dis- 
tribution shown in figure 5 (left panel) of Purcell et al. (2007). Grey hatching 
and the black point with errorbars indicate likely values for the Milky Way 
(Smith et al. 2007; Li & White 2008; Bell et al. 2008; McMillan 201 1) and 
M31 Watkins et al. (2010); Courteau et al. (2011) respectively. 



ciency per unit halo mass increases steeply along the halo mass 
function, up to a maximum around the Milky Way mass (e.g. 
Moster et al. 2010). Because even the most massive accreted haloes 
typically have much lower galaxy formation efficiency than the 
main halo, A/ acc remains a small fraction of M* and depends 
strongly on A^oop 2 ] For M200 > 10A/ pca i;, the most massive pro- 
genitor haloes have a galaxy formation efficiency comparable to or 
even higher than the main halo and Af acc makes up a much larger 
fraction of M+. However, as the M* -A/200 relation flattens, the 
number of progenitors continues to scale with A/200, but the stel- 
lar mass per massive progenitor does not. The result is that M acc 
scales more slowly with A/200 than it does below A/ pca k. 

The Gil model follows the entire hierarchy of galaxy forma- 
tion and with particle tagging we have now computed the dynami- 
cal evolution of all satellite disruption events in that model directly 
from Millennium II. These techniques mean that our model is sig- 
nificantly more accurate and detailed than the empirical 2 = scal- 
ing relations used by Purcell et al. (2007) for A/* (A/200 ) , R50 (M t ) 
and the disruption of satellite galaxies, so it is useful to revisit their 
analysis. Fig.[l5]presents the data from Fig.[l4]in the same way as 
figure 4 of Purcell et al. (2007), showing the ratio of M acc to 
as a function of A/200 • As expected, our results are qualitatively 
similaf^jto those of Purcell et al.. However, our model predicts a 



12 In the extreme case, sharp thresholds for galaxy formation at very low 
M200 (the atomic hydrogen cooling limit, reionization etc.) will result in 
central galaxies that do not accrete any luminous progenitors and have 
A/ acc = 0. 

13 Purcell et al. (2007) define their virial quantities at an overdensity of 



significantly steeper relation. An average accreted stellar mass frac- 
tion of 30 per cent is reached in haloes with A/200 ~ 1O 12 ' 5 M 
rather than A/200 ~ 1O 13,5 M . 

Comparing these results to observations is not straightforward, 
because of the difficulty of separating accreted and in situ stars. In 
both our results and those of Purcell et al. A/ acc does not represent 
either the stellar halo mass or the total spheroid mass (halo plus 
bulge). This is because an unknown fraction of spheroid stars (both 
bulge and halo) may form in situ. In the Milky Way, a rough lower 
limit for A/ acc is the total mass of halo stars (e.g. Bell et al. 2008), 
approximately 1 - 2.5 per cent of the mass of the disk and bulge 
(e.g. McMillan 201 1). Estimates of A/200 for the MW using differ- 
ent methods range over an order of magnitude, from ~ 4 x 1O 11 M0 
(18th percentile of Smith et al. 2007) to ~ 4 x 1O 12 M (the 95th 
percentile of). 

Given these ranges, and making the naive assumption that all 
halo stars are accreted while all disk and bulge stars formed in situ, 
the grey shaded region in Fig. [15] illustrates the region of likely 
Macc/M* for the Milky Way. For A/200 > 1O 12 M this region 
lies below the median curves of the two models, which both pre- 
dict that ~ 3 per cent of the total stellar mass in a dark halo of mass 
A/200 ~ 1O 12 M has been accreted. This hints that the Milky Way 
has a less massive stellar halo than average, or else that the lowest 
estimates of A/200 are preferred. It could also indicate that A/ acc for 
the MW should include substantially more mass than the observed 
stellar halo, which seems plausible as the density profile of the stel- 
lar halo inside the Solar radius and its relationship to the bulge and 
thick disk are unknown. On the other hand, the Milky Way could 
soon substantially increase the mass of its stellar halo by disrupting 
the LMC. Such stochasticity in the time of accretion of the most 
significant halo progenitor may be an important contributor to the 
scatter in Fig. 1 1 5 1 

Similar data are available for M31; its stellar halo mass frac- 
tion has been estimated by Courteau et al. (2011) and M200 has 
been constrained to be similar to that of the Milky Way by the 
timing argument ( ) and satellite kinematics (Watkins et al. 2010). 
These data are represented by a black point with error bars (2cr in 
A/200 from Watkins et al. 2010 and ±50 per cent in accreted mass 
fraction from Courteau et al. 2011). This point agrees well with 
both our model and that of Purcell et al. (2007). The accreted mass 
fraction depends on the assumed density profile for the stellar halo 
of M31, which in turn depends on the bulge/disk/halo decomposi- 
tion. Because a substantial amount of light in the spheroid of M3 1 
is known to be associated with a single tidal stream, its measured 
mass may be closer to the definition of A/ acc than the Milky Way 
stellar halo. 

Both models therefore appear roughly consistent with the 
(limited) data on late type galaxies in haloes of A/200 ~ 1O 12 M . 
They disagree with one another at higher mass, where most galax- 
ies are ellipticals and even a lower limit to the accreted component 
is hard to identify by decomposing the light profile. Our model pre- 
dicts that a late type galaxy in a halo of 1O 13 M should have ~ 30 
per cent of its total stellar mass in an extended n ~ 3-4 spheroid 
of accreted stars (Purcell et al. predict ~ 10 per cent). This can be 
tested with deep images of nearby massive disk galaxies (including 
SOs) for which the total mass is well constrained by the rotation 
curve. 



A = 337 rather than A = 200 as we do. We have not corrected their 
results for this difference in Fig. |15| doing so would shift their curves to the 
right (higher virial mass) by < 0.1 dex for plausible NFW concentrations. 
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Finally, we discuss the progenitor galaxies that contribute stars 
to the accreted component in our models. A progenitor is defined as 
a galaxy that is disrupted within the 'main branch' of a halo merger 
tree - each progenitor may have many progenitors of its own, but 
all these are grouped together in this definition based on only the 
final level in the hierarchy. The left-hand column of Fig. |16| shows 
the mass of the most massive progenitor in each of our galaxies. 
The upper and lower panels of this plot separate two classes of 
progenitor, according to where their stars settle in the main halo 
after they have been accreted. We define 'bulge' progenitors to be 
those that deposit more than half their stars within a radius of 3 kpc 
from the centre of the main halo at z — 0. The rest are classified as 
'stellar halo' progenitors (the same definition was used in CIO). 

The most massive 'bulge' progenitor is typically the most 
massive of all accreted galaxies, which is not surprising because 
more massive satellites suffer more dynamical friction and sink 
quickly to the centre of the potential. Thus the relation between 
the mass of the most massive bulge progenitor and the mass of the 
dark halo is similar to that shown in Fig. [14] with a steep slope for 
late-type galaxies and an approximately constant value for ellipti- 
cals. On the other hand, the mass of the most massive stellar halo 
progenitor shows the same trend with M200 in both late types and 
ellipticals, thus our model predicts the outer stellar haloes of all 
galaxies to be equally diverse at a fixed M^oo- The most massive 
bulge progenitor is ~ 10 times as the most significant contribu- 
tor to the halo up to M200 ~ 10 13 Mq; at higher M200 the most 
significant bulge and stellar halo progenitors have similar mass. 

The right-hand column of Fig.[l6]shows the mass of the most 
massive bulge and halo progenitors as a function of the total ac- 
creted stellar mass of the main galaxy in the same region. The ac- 
creted bulges of late type galaxies usually acquire at least ~ 40 
per cent of their stars from one progenitor, reaching > 90 per cent 
in many cases. On the other hand, the most massive contributor 
to the accreted stellar halo typically accounts for no more than 
~ 30 per cent of its total mass, and rarely exceeds 50 per cent. In 
more massive galaxies (M* > 10 11,2 Mq, which mostly have Gil 
B/T > 0.9) we find very few cases where more than 70 per cent 
of accreted bulge or halo stars originate in a single progenitor. 

We can quantify the number of significant progenitors of the 
accreted component as a whole with the definition (see C10) 



N sig = 



(2) 



where m p is the stellar mass of progenitor p. This expression is 
equal to unity when one progenitor contributes all of the accreted 
mass, and N when N progenitors contribute equal mass. As dis- 
cussed above, progenitors with M200 > M poa k (the most efficient 
halo mass for galaxy formation) have roughly the same stellar mass 
as the central galaxy they are accreted on to. This means that mas- 
sive ellipticals can undergo mergers that are 'minor' in terms of 
M200 (and thus numerous) but 'major' in terms of accreted stellar 
mass. We find the most massive galaxies in our sample typically 
have ~ 10 equally important progenitors to their accreted com- 
ponent (see Laporte et al. 2012). For ellipticals there is a reason- 
ably tight correlation between the number of significant progeni- 
tors and Af2oo; the lowest mass ellipticals in our sample typically 
result from a single major merger and thus have only one signifi- 
cant progenitor which contributes a large fraction of its mass to the 
accreted bulge. Late type galaxies have a larger scatter in the num- 
ber of significant progenitors of A/ acc . They can be dominated by 
one massive object, in agreement with the findings of C10, but we 
also find cases with 10 or more significant progenitors. 
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Figure 16. The stellar mass of the most massive progenitor, M*. mm , as 
a function of total halo mass M200 (left) and the fraction of stellar mass 
acrreted from the most massive progenitor as a function of galaxy stellar 
mass JVf* (right). Upper and lower panels separate progenitors by the half 
mass radius r^o of their stellar debris at 2 = 0: rgn < 3 kpc (the 'bulge') 
and rso > 3 kpc (the 'stellar halo'), respectively. The dashed black lines in 
the left-hand panels show the mass corresponding to 1 X 10 — 4 A/200- 



8 CONCLUSIONS 

We have used dark matter particles in the Millennium II simulation 
as dynamical tracers of stellar populations, in order to study the 
hierarchical assembly of galactic structure in the CDM model. Our 
most important conclusions can be summarised as follows: 

(i) The stellar mass surface density profiles of galaxies in our 
model, E(-R), are well described by the sum of two Sersic models, 
corresponding to the separate contributions of stars formed in situ 
and stars accreted from other galaxies. 

(ii) The surface density of in situ stars falls off more rapidly 
with radius than the accreted component. In situ stars only make 
a significant contribution to surface brightness profiles out to R ~ 
10 kpc in M200 ~ 1O 12 M halos or 40 kpc in Af 2 oo ~ 1O 14 M 
haloes for / m b ~ 0.01. 

(iii) The outer isophotes of all massive galaxies are dominated 
by accreted stars, which extend to the virial radius in most systems. 

(iv) E(.R) shows very little scatter from galaxy to galaxy at fixed 
M200, particularly in cluster-scale haloes (M200 > 10 13 ' 5 Mq). 
This is the consequence of strong correlations between halo mass, 
central star formation efficiency and the galaxy progenitor mass 
function, and is a basic feature of galaxy formation in the CDM 
model (e.g. Purcell et al. 2007). Comparison at fixed M* shows 
more scatter, because of scatter in the relationship between M* and 
Af 2 oo- 

(v) E(ii) profiles can be in situ dominated or accretion dom- 
inated. Accretion-dominated profiles are more extended and have 
higher Sersic index compared to in situ-dominated profiles. They 
have approximately 'power law' profiles from 10-100 kpc that 
show no clear inflection at the transition between in situ and ac- 
creted components and are associated with haloes that host ellip- 
tical galaxies, according to the criteria described by Gil. These 
massive haloes are usually subject to 'major' mergers (violent re- 
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laxation) after z ~ 1, which reshape their in situ component. They 
also suffer strong suppression of star formation by AGN feedback 
which prevents a compact core of in situ stars forming at low red- 
shift. 

(vi) In situ dominated profiles are typical of less massive haloes 
(up to M200 ~ 10 13 Mq). In situ stars have a compact exponen- 
tial distribution by construction in our model. Extended (R 5 o ~ 
10 kpc) accreted spheroids begin to dominate only at R > 10 kpc, 
causing a clear inflection in the circularly- averaged surface bright- 
ness profile. Neither compact bulges nor compact exponential el- 
liptical galaxies are created by accretion in our model. 

(vii) The transition from in situ dominated profiles to accretion 
dominated profiles with increasing M200 is less clear when studied 
in terms of central stellar mass M+, because of scatter in the M*- 
M200 relation in our model. 

(viii) Sersic parameters describing the surface density of ac- 
creted stars have roughly power-law scalings with M200. At fixed 
M200, the surface density profiles of accreted stars have larger half- 
mass radii i?5o,acc (~ 10 kpc) and Sersic index n acc (~ 3-4) than 
those of in situ stars. 

(ix) For galaxies with B/T > 0.9 in our semi-analytic model, 
the accreted component scales approximately as i?so,acc oc M^. 
For galaxies with B/T ^ 0.9, -Rso.acc scales more slowly with 
M200, such that for M200 ~ 10 12 Mq, the accreted components 
of elliptical galaxies are smaller and more concentrated than the 
stellar haloes of late type galaxies. (When the total light profile is 
considered, late type galaxies in our low-mass dark haloes are nev- 
ertheless smaller than our 'ellipticals', because they are dominated 
by a compact core of in situ stars.) 

(x) Intracluster light can be explained through tidal stripping, 
analogous to the stellar haloes of Milky Way-like galaxies. We find 
no evidence for a separate ICL component in our models; rather 
the ICL is the outward extension of accreted stars that make up 
the dominant fraction of the central light. However, this does not 
mean that ICL stars have the same progenitors (hence chemistry 
and kinematics) as accreted stars in the centres of galaxies. 

(xi) Radial variations in the stellar population of the accreted 
component can occur in cases where different progenitors domi- 
nate the light profile at different radii. In our model, the accreted 
spheroids of massive cluster central galaxies typically have ~ 10 
significant progenitors with similar masses and similar (old) stellar 
populations. This results in weak population gradients for the ac- 
creted component as a whole. (The two-dimensional light distribu- 
tion could show other, perhaps stronger, differences - for example 
in the radial variation of isophote ellipticity and position angle.) 

(xii) We have stacked SDSS images of galaxies in bins of stellar 
mass to obtain average surface density profiles that can be com- 
pared directly to the results of our simulations. We find a weak 
trend of surface density amplitude with stellar mass but no clear 
change in profile shape. This is consistent with our simulations 
when we consider biases in the measurement of total stellar mass. 

(xiii) We have also compared the results of our model with data 
from the literature, including deep surface brightness profiles of in- 
dividual galaxies and stacks of LRGs and BCGs. We find qualita- 
tive agreement the shape and amplitude of the stellar mass surface 
density profile and its variation from low mass to high mass dark 
matter haloes. 

The methods we have presented in this paper are approximate. 
We have identified several important areas of disagreement with the 
data and limitations that could be improved on in future work: 

(i) Our free parameter / m b is only weakly constrained by a com- 



parison to the observed galaxy size-mass relation, particularly at 
high stellar masses where the observations are most affected by 
systematics in the measurement of faint isophotes. A similar con- 
straint comes from our comparison to stacked surface density pro- 
files from SDSS; in both cases, we find acceptable values in the 
range 1% < / m b < 5%. An improvement on the method would 
be to remove the freedom of choosing / mb by computing it on the 
basis of other physical quantities in our semi-analytic model. In 
this case / m b would no longer be constant. However, such com- 
plications may not be justified given the severe approximation of 
neglecting the gravity of stars and cold gas. 

(ii) Our model for the distribution of in situ stars is a very crude 
representation of the inner regions of galaxies. Apart from be- 
ing unable to represent disks, we cannot properly resolve compact 
bulges formed by dissipation. Furthermore, bars, pseudobulge and 
lens structures are known to exist in real galaxies and are thought to 
be the result of secular processes that are not included in our sim- 
ulations. These aspects of galactic structure are much better mod- 
elled by high resolution hydrodynamic simulations. Such simula- 
tions could be used to test our highly approximate model for the 
dynamics of in situ star formation. 

(iii) There is evidence in real data for an excess of light over our 
model in massive galaxies at > 100 kpc, particularly for cluster 
central galaxies (e.g. Zibetti et al. 2005). The GIMIC SPH simu- 
lations (Crain et al. 2009; Font et al. 2011) show shallower halo 
slopes at large radii for Milky Way-mass galaxies. As all stars at 
these radii are accreted, this suggests that our model is deficient ei- 
ther in the treatment of the orbits of galaxies that are being stripped 
or the treatment of the stripping process itself. For example, it may 
be the case that the angular momentum of disk galaxies affects the 
orbital energy distribution of their tidal tails, or that ram-pressure 
stripping of cold gas leads to more rapid disruption of satellites; 
in both cases, more stripped stars would be deposited on weakly 
bound orbits with large apocentres. 

The next generation of deep sky surveys (culminating in 
LSST) will reach the low surface brightness limits required to de- 
tect stellar haloes and accretion remnants around the majority of 
low-redshift galaxies. The dawn of this new era in the observa- 
tion of galaxy structure is a strong motivation for the further study 
of these regions in simulations. We have shown that the particle tag- 
ging technique is a straightforward and well-constrained extension 
of the semi-analytic method with a number of interesting appli- 
cations that merit further investigation. We hope to explore some 
of these in future work, including stellar population gradients in 
early-type galaxies; the kinematics of diffuse light; the effects of 
interactions on the structure of satellite galaxies and environmental 
trends; the frequency of tidal features and their correlations with 
other galaxy properties; and the intracluster light of massive clus- 
ters. 
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APPENDIX A: DIFFERENCES WITH CIO AND 
NUMERICAL CONVERGENCE 

There are a number of minor differences between the tagging pro- 
cedure we use on Millennium II and that used on the Aquarius sim- 
ulations by CIO. First, as discussed at length in the main text, we 
now tag particles to represent in situ stars in our 1872 target cen- 
tral haloes as well as those formed in their progenitors and satellite 
galaxies. In CIO, the distribution of in situ stars in the main halo 
was not considered, but only because they were not the focus of 
that paper and excluding them made the computation less demand- 
ing. 

Second, the semi-analytic component of our model must deal 
with galaxies in dark matter haloes that fall below the mass reso- 
lution limit of the simulation because of tidal stripping. We do this 
in the same way as Gil, by estimating the time for inspiral of the 
satellite due to dynamical friction in its parent halo, and merging 
it with the central galaxy after that time. The satellite can also be 
disrupted sooner if its density falls below that of the dark halo at 
the satellite's pericentre. This approach differs from CIO, where 
galaxies in unresolved subhalos were merged instantly with the 
main galaxy in their parent halo (in which case the disruption of 
the N-body subhalo is equivalent to the disruption of the galaxy, 
as it would be in an SPH simulation with collisionless star parti- 
cles). CIO did this because of the high resolution of their Aquar- 
ius simulations, which meant that it made almost no difference to 
their results. However, at the lower resolution of Millennium II, 
allowing galaxies to survive after the disruption of their N-Body 
subhalos is necessary for convergence of the semi-analytic model 
and its agreement with observations (details are given in Gil; see 
also Font et al. 201 1). Therefore, our tagging method must make al- 
lowance for new stars forming in galaxies with unresolved haloes. 
We do this by tagging those stars to a single particle, the most- 
bound particle of the halo at the time it was last resolved. Because 
very few stars form in such haloes, this makes no practical differ- 
ence. 

There is a more significant issue related to the treatment of 
semi-analytic galaxies with unresolved haloes. Our analysis (for 
example, in constructing density profiles) treats the stars in these 
galaxies as having been stripped, because their dark matter par- 
ticles are bound to the main dark matter halo. This includes the 
most-bound particle to which we tag any residual star formation. 
Again, this treatment is similar to conventional hydrodynamic sim- 
ulations, which usually do not track galaxies below the resolution 
limit. However, because of this, our definitions of stellar mass are 
not consistent between the semi-analytic model (where these all 
stars belong to a galaxy) and our tagged-particle model (where they 
are all stripped). These ambiguous stars are easy to identify in the 
model, and we have verified that this choice makes no significant 
difference to our results. 

A related problem is that the lower resolution of Millennium II 
(relative to the Aquarius simulations of CIO) makes all satellite 
galaxies somewhat easier to disrupt (at fixed mass), because their 
density is artificially reduced at radii below the force softening 
scale. Finally, it may be the case that Millennium II does not re- 
solve a significant number of faint progenitors at all, particularly 
for the least massive main galaxies in our sample. However, as the 
much higher resolution simulations of CIO have shown, most of 
the stellar mass accreted by Milky Way mass haloes comes from a 
small number of their most massive progenitors, which are well re- 
solved in Millennium II. In Fig. |Al| we compare the high-resolution 
Aquarius simulations of CIO with the same haloes simulated at 
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the resolution of Millennium II. We find that resulting stellar mass 
density profiles have converged above the softening scale of Mil- 
lennium II. The details of profile shape that distinguish different 
Aquarius haloes in C 10 are reproduced in Millennium II, which has 
a particle mass almost four orders of magnitude larger. This is the 
case for both the accreted and in situ components. This implies that 
the resolution effects described above are not very important for 
our results (and notably, that the underlying semi-analytic model of 
Gil has also converged with regard to the predicted stellar mass of 
the Aquarius progenitors). 

The final difference with CIO is that we do not postpone the 
tagging process to later snapshots in cases where the target halo 
is deemed to be out of equilibrium. Although well-motivated, CIO 
found that this procedure makes the implementation of the method 
much more complex but has little influence on the outcomj^] 
Other technical subtleties, including tagging a fixed number of par- 
ticles rather than a fixed fraction in subhalos that are losing mass, 
are dealt with as described in CIO. 



APPENDIX B: NUMERICAL APPROXIMATION FOR 
IN-SITU SIZES 

In Section [23] we used a simple numerical approximation to our 
particle tagging results to show how the / m b parameter is related 
to the scale of in situ stars (at least in haloes that are not per- 
turbed by mergers after their main epoch of star formation). This 
approximation was based on computing the density profile of mass 
more bound than a given energy in an idealized NFW potential, for 
which we used the following straightforward computational proce- 
dure. We used scipy 0.9.0 to solve the necessary numerical 
integrals (with scipy . integrate . quad). 

Appendix A of Widrow (2000) gives an analytic fitting for- 
mula for the distribution function of an isotropic NFW halo, F(e), 
where e oc — (E — $oo) is the normalized energy and ^ oc 
— ($ — $oo) the normalized potential. In these expressions, E — 
v 2 jl + 4>(r) is the total energy in physical units (v is velocity); 
$(r) is the Newtonian NFW potential, and "Soo is the limit of $(r) 
as r — > oo. In the convention used by Widrow (2000), which we 
follow, these energies are normalized such that ^(0) = 1. The di- 
mensionless formula for F(e) can be scaled to any halo given the 
parameters of its NFW profile (e.g. virial mass and concentration). 
These parameters are measured in Millennium II using the SUB- 
FIND algorithm (specifically, we measure M200 and i? max , the ra- 
dius at which V c 2 (r) = GM(< r)/r peaks; for this approximate 
method we assume all haloes are described by a spherical NFW 
profile). 

Given F(e) from the Widrow (2000) formula, the differen- 
tial energy distribution is given by dM(e)/de = F(e)g(e), where 
g(e) is the 'density of states' at a particular energy (equation 1 1 of 
Widrow 2000). Evaluating g(e) as an integral over radius requires 
the radius of maximum excursion for a given e, r max , such that 
^(w(f)) = 0. For a range of e, we find this root numerically 
(with scipy . optimize . brent q) using an analytic expression 
for W(r) (e.g. table 1 of Widrow 2000) and then evaluate g(e) by 
numerical integration. 

With these expressions for F{e) and g(e), we tabulate val- 
ues of dM/de over the range e — [0.1, 1] and fit them with a 

14 At the level of the other approximations in this method, it is arguable 
that allowing assignments to non-equilibrium haloes might be a reasonable 
representation of the messy nature of star formation in mergers. 



cubic spline. We integrate this spline to tabulate M(> e), the 
cumulative mass more bound than e. From M(> e) we deter- 
mine e m b, the cut-off value of e corresponding to a particular 
/mb = M(> e m b) /A/200 • The final step is to obtain the density 
profile of all mass with energy less than e m b by evaluating equa- 
tion 4 of Widrow (2000): p(r) oc /* F(e)(* - e) 1/2 de. 



APPENDIX C: SDSS STACKING ANALYSIS 

In Section|6]we stacked galaxies from SDSS DR9 for comparison 
with our models. This appendix describes our method for selecting 
galaxies in SDSS and stacking their images. We intend to explore 
a number of important systematic uncertainties in more detail and 
present further results in a separate paper (D'Souza et al. in prepa- 
ration), so our analysis here should be considered preliminary. 

Our starting point is the MPA-JHU SDSS 'value-added' cata- 
logue, which provides an estimate of stellar mass for galaxies with 
spectra in DR7 based on fitting an SED to their modelMag pho- 
tometry (Kauffmann et al. 2003a; Salim et al. 2007). From this cat- 
alogue, we selected isolated central galaxies in the redshift range 
0.07 < z < 0.08 by applying the criteria of Wang & White (2012): 
a galaxy of apparent r band magnitude m is considered isolated if 
there are no galaxies in the spectroscopic catalogue at a projected 
radius 7? < 0.5 Mpc and velocity offset \5z\ < lOOOkms -1 
with magnitude m' < m + 1, and none within R < 1 Mpc and 
\5z\ < lOOOkms" 1 with m' < m. We make no selection on 
colour or morphological type. 

We created 1 Mpc 2 mosaics in the g, r and i bands centered 
on each galaxy in our sample using the 'corrected' sky-subtracted 
frames from the SDSS Data Release 9 image server and SWarp 
(Bertin et al. 2002). These three mosaics were stacked together 
to make a 'master image' from which a 'master mask' was ob- 
tained using SExtractor (Bertin & Arnouts 1996). Other galax- 
ies in the field were conservatively masked by convolving the 
master image with an 8 x 8 pixel top hat kernel before running 
SExtractor. The master mask was applied to each individual 
mosaic. The masked mosaics were then transformed to z — 0.08 
with the flux-conserving IRAF task geotran, cropped to a uni- 
form size of 1200 x 1200 pixels (690 x 690 kpc at z ~ 0.08) and 
corrected for extinction following Schlegel, Finkbeiner & Davis 
(1998). 

We assume that the sky subtraction provided by the SDSS 
imaging pipeline is adequate for our analysis. The pipeline sky sub- 
traction was known to have shortcomings in early data releases, but 
was revised in DRfJ^Jto improve the photometry of extended low- 
surface brightness regions around low redshift galaxies (Blanton 
et al. 201 1). We confirmed the quality of the sky subtraction in the 
DR9 images by carrying out tests using the earlier DR7 images 
with our own background subtraction, and also by comparing DR9 
surface brightness profiles for galaxies in the Virgo cluster with the 
results of Kormendy et al. (2009). 

We binned our SDSS sample for stacking in mass bins as we 
did our simulations in Fig. [T2] using the mode of the MPA-JHU 
mass PDF corrected to a Hubble parameter h = 0.73 for consis- 
tency with the cosmology of our simulations"*] This resulted in 



See 



www . sdS53 . org/ cir 9 /imaging/ images .php| and |www . | 



|sdss3 . org/dr9/imaging/caveats .php| 
±u We have not accounted for a possible overestimate of the DR7 
modelMag magnitudes owing to uncertainties in the background subtrac- 
tion of extended galaxies pre-DR8. 
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~ 600-1200 SDSS galaxies per bin. Mosaics in the gri bands 
were stacked separately using IRAF imcombine, taking the mean 
valu^Jof each pixel after clipping at the 10th and 90th percentiles. 
Images were centered before stacking but were not axis-aligned. 
We also removed a small residual background from each stack such 
that the stellar mass surface density falls to zero at the periphery of 
the mosaic. The g, r and i band mosaics for each galaxy were not 
included in the corresponding stacks if they included stars brighter 
than r = 12, if they fell in the upper 10 th percentile of the dis- 
tribution of pixel RMS for each mosaic in their mass bin, if they 
had more than 75% of their central 60 arcseconds masked, or if the 
masking algorithm failed due to a crowded field. 

In a given mass bin, the azimuthally averaged mass profile 
was derived from the r band stack assuming a constant mass-to- 
light ratio. The M/ L r assumed for each bin is the average of M/L 
for the galaxies in the bin, derived from their MPA-JHU mass and 
r band modelMag; in order of M* for the four bins, M / L r — 
[2.154, 2.240, 2.391, 2.516]. We also used the g, r and i stacks to 
derive a radially varying M/L r profile based on the colour of the 
light in each annulus, using the prescriptions of Bell et al. (2003) 
with a Chabrier IMF. This only makes a significant difference in 
the outermost radial bins of each profile. However, the potential for 
large colour errors at faint magnitudes (particularly those involving 
the i band) adds considerably to the uncertainty in these radially 
varying profiles. 

We have not carried out deconvolution of the PSF at any stage 
of our analysis. Galaxies in the mass and redshift range of our sam- 
ple are typically well resolved, so the PSF mostly affects the very 
inner part of the light profile. We have compared the profiles of 
individual galaxies in our stacks with similar galaxies at lower red- 
shift, with SDSS profMean profiles and with deep data from the 
Virgo cluster (Kormendy et al. 2009) to verify that the PSF only 
has a significant effect (very roughly, of the order of > 0.1 mag- 
nitudes per square arcsecond in the surface brightness profile) at 
R < 5 kpc. Nevertheless, particularly for the i band magnitude 
and the measurement of colours, PSF effects are known to be sig- 
nificant at much larger radii (de Jong 2008) and a complete analysis 
should include a more thorough quantification of the PSF (e.g. Tal 
& vanDokkum 2011). 

We estimated the uncertainty in each annulus in the constant 
M/L case as the sum in quadrature of Poisson error in the flux 
per contributing pixel and the average RMS pixel-by-pixel devia- 
tion of each image from the stack in that annulus, the latter term 
accounting for the sample variance. A bootstrap estimate of vari- 
ance would be preferable to understand the effects of sample vari- 
ance, but proved to be computationally expensive. Bootstrapping 
on sub-samples of 100 galaxies suggested that our algorithm un- 
derestimates the variance in the outer regions of the profile by a 
factor of 2, so we multiply our variance estimate by 2 to obtain the 
error bar in Fig[T3] This crude estimate of uncertainty is sufficient 
to indicate the largest radius to which each stacked profile is robust, 
but it could be made substantially more accurate with further work. 

In our selection of SDSS galaxies for stacking, we have used 
the MPA-JHU stellar mass estimates. These masses are derived 
from the galaxy modelMag magnitude, which attempts to correct 
for undetected light in regions of low surface brightness by fitting 
the observed surface brightness profile to one of two analytic mod- 
els (exponential and r 1 ^ 4 ) and adopting the total magnitude of the 
best fit. However, because the fits are to truncated light profiles 



Median and mode stacking produced very similar results. 



(owing to the finite surface brightness limit, of more importance 
for more extended galaxies) and assume a fixed Sersic index n = 4 
for the early-type model, they are likely to suffer from a bias rela- 
tive to the true total light similar to that of the Petrosian magnitude 
petroMag. In theory petroMag is a substantial underestimate 
of the total light for galaxies with Sersic index n > 4, but includes 
almost all of the light for galaxies with n ~ 1 (Graham et al. 2005; 
Lauer et al. 2007; Blanton et al. 201 1). 

The bias in modelMag is more difficult to reproduce in our 
simulations because it depends on the fitting process itself (for ex- 
ample, in the behaviour of fits to multi-component galaxies, or with 
substantial noise). Furthermore, we find that the median offset be- 
tween modelMag and petroMag for galaxies in our SDSS sam- 
ple is generally small (~ 0.1 magnitudes). We therefore compute 
the 'Petrosian' mass from the simulations and use this as a proxy for 
the SDSS modelMag mass in Fig.[T3] We do so with an algorithm 
analogous to that used for Petrosian flux reported by SDSS (i.e. 
the mass within 2r p , using the same definition of r p ). The result- 
ing Mp C t greatly underestimates Ah for the most massive galaxies, 
thereby redistributing galaxies with high-n profiles to lower mass 
bins in our simulated stacks. 
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